Response to Perturbation During Quiet Standing Matches Delayed State Feedback With Precisely Tuned Control Gains


 Postural sway is a result of a complex action-reaction feedback mechanism generated by the interplay between the environment, the sensory perception, the neural system and the musculation. Postural oscillations are complex, possibly even chaotic. Therefore fitting deterministic models on measured time signals is ambiguous. Here we analyse the response to large enough perturbations during quiet standing such that the resulting responses can clearly be distinguished from the local postural sway. Measurements show that typical responses very closely resemble those of a critically damped oscillator. The recovery dynamics is modelled by an inverted pendulum subject to delayed state feedback and is described in the space of the control parameters. We hypothesize that the control gains are tuned such that (H1) the response is at the border of oscillatory and nonoscillatory motion; (H2) the response is the fastest possible; (H3) the response is a result of a combined optimization of fast response and robustness to sensory perturbations. Parameter fitting shows that H1 and H3 are accepted while H2 is rejected. Thus, the responses of human postural balance to “large” perturbations matches a delayed feedback mechanism that is optimized for a combination of performance and robustness.


Introduction
For over 100 years, responses to perturbations have been used to investigate the feedback control of human balance [1][2][3][4] . These studies established that human balance is not maintained by stereotyped reflexes. Instead, with development, balance control emerges as the nervous system learns to apply generalized rules for maintaining balance 5 . In healthy individuals this ability to adapt and improve balance in a feedback-driven manner does not appear to decline with age 6 . This observation underscores the current development of perturbation-based balance training protocols to improve reactive balance in the elderly as a way to reduce their risk of falling 4,7 .
A large variety of perturbations have been used to disturb human standing balance including sudden platform translations, pulls and tugs 2, 3, 8-11 . The nervous system responds with a continuum of "ankle" and "hip" strategies, the exact combination depending on the trade-offs between the required effort and the degree of postural instability to be overcome [12][13][14][15] . Despite these observations, theoretically-motivated investigations suggest that common underlying principles may be at work 9,[15][16][17][18] . Thus, the challenge has become to identify the nature of the governing principles 8,19,20 .
The dynamics of human postural sway during quiet standing with eyes closed is very complex and has been described in terms of stochastic and even chaotic motions 21,22 . The underlying mathematical model can be either an inherently stochastic process 23,24 or a deterministic nonlinear feedback mechanism governed by some kind of intermittent control [25][26][27][28] or their combinations [29][30][31] . One of the simplest physiological mechanisms for generating a chaotic motion is the interplay between a time-delayed, sampled data system and a sensory dead zone, i.e., corrective actions take place only when the sensory inputs exceed some threshold values 22 . Model based analysis of postural sway is therefore a difficult task since the governing deterministic dynamics may be hidden in the seemingly noisy/chaotic response. Here we employ perturbations in the anteroposterior (AP) and the posteroanterior (PA) directions during quiet standing that are large enough to produce excursions that are significantly larger than the magnitude of the fluctuations in postural sway and the size of sensory dead zones. In this way, the response to the perturbation is not affected by the local noisy/chaotic dynamics and the underlying feedback mechanisms can be identified distinctly by numerical fitting techniques.
Inverted pendulum models are widely used to investigate human balance 20,[32][33][34][35][36] and are currently thought to be "functionally correct" 35,37 . Reactive balance control, that is maintaining balance in the response to a perturbation, is inherently a feedback sensorimotor process in which muscles are activated in direct response to task-level error 15,[38][39][40] . However, the time delayed nature of the feedback control has important implications for the feedback controlled response of balance to perturbations. Our concept is illustrated in Figure 1. This figure shows the stability diagram for an inverted pendulum stabilized by a time-delayed proportional-derivative (PD) feedback controller. The stable region in the plane (p, d) of the control gains is characteristically D-shaped 32,41 . For every choice of the control gains located within this D-shaped region, the upright position recovers from a perturbation; however, not all control gains exhibit the same dynamical response to the perturbation. The perturbed responses range from a monotonic exponential recovery to the upright position to recoveries which exhibit an oscillatory component. The boundary between the monotonic and oscillatory responses is formed by the node-spiral separation line (black-green dashed line). The exponential decay of the response is shown by solid contour lines: the smaller the value of γ 1 , the faster the response. Namely, |θ (t)| θ 0 e γ 1 t where θ 0 is the initial angle. γ 1 < 0 is called exponential decay rate 42 . The fastest responses to perturbations are for choices of the gains in the lower left quadrant of the stability diagram. The optimal point with respect to the response's settling time is (p * , d * ) (black × marker), which is associated with the maximal achievable decay rate γ 1 = γ * . This parameter point lies on the node-spiral separation line and divides it into two sections indicated by black and green color. Small changes in p or in d in the lower (black) branch of the separation line results in significantly larger changes in the decay rate γ 1 than the same changes do in the upper (green) branch. Hence, the system is more robust to changes in the control gains when operating at a parameter point on the upper (green) branch rather than on the lower (black) one. Note that perturbation of the gains can directly be linked to perturbations in sensory perception since the actual control force is determined as the product of the control gains p and d and the corresponding sensory inputs, the angle θ and the angular velocityθ . This suggest that the control gains are rather to be selected based on a goal function which takes into account both fast response and robustness. This would give an operation point (p,d ) on the upper branch indicated by green × marker.
The main goal of this paper is to compare the dynamics of the time-delayed PD feedback model to measured responses to perturbations during quiet standing. Control gains and feedback delays are estimated by fitting the response of the mechanical model to the measured time histories. The fitted parameter point is indicated by (p,d) in Figure 1b (black marker). We pose three hypotheses related to the location of the fitted control gains.

H1
The fitted control gains are tuned towards the node-spiral separation line, which would indicate that reducing oscillatory response is one of the main goal of the feedback mechanism.
H2 The fitted control gains are close to the parameter point (p * , d * ) associated with the fastest possible recovery of balance to perturbations. H3 The fitted control gains are located close to a parameter point (p,d) on the upper (green) branch of the node-spiral separation line (i.e.,p > p * andd > d * ) that assures fast recovery even upon small uncertainties in the control gains.
We base our hypotheses on three observations: 1) PD controllers with gains located in the lower left quadrant of the stability diagram are most robust to the effects of random perturbations 43 ; 2) expert pole balancers increase maneuverability while minimizing energetic costs for balance control by adapting gains in the lower left quadrant of the stability diagram 44 ; and 3) the control gains for subjects who do ball-and-beam balancing are progressively tuned towards the node-spiral separation line as their skill improves 45 . Statistical analysis show that H1 and H3 are accepted and H2 is rejected.

Experiments and responses to perturbation
The measurement setup, a typical response and measured responses for AP and PA perturbations are shown in Figure 2. After the abrupt release of the weight at time instant t 0 , the tilt angle θ first reaches a maximum then recovers to the angle that corresponds to normal posture. The oscillations before the perturbation (t < t 0 ) and after recovery (t > t 2 ) were significantly smaller (std 0.15 • in average) than the maximum excursion caused by the perturbation (3.87 • in average). Since the force was released at an unexpected instant 46 , the subject corrective action was delayed by the reaction time τ, which, in theory, is the difference between the initial time t 0 and the time instant where the response has inflection point, i.e., the angular acceleration θ changes sign. Responses were typically free of overshoot and resemble the responses of a critically damped oscillator 47,48 .
Recorded time signals were used to fit parameters p, d and τ in the mechanical model.

Identified control parameters
The control gains p and d and the time delay τ were estimated for all the 20 trials (10 PA and 10 AP) for all the 10 subjects. The stability charts together with the identified control gain parameters are shown in Figure 3 for each subject. The stable region (thick black solid curve) and the node-spiral separation line (dashed black-green curve) correspond to the average delay of the overall 20 trials per subject (the average delay is indicated in each panel). The larger the average delay, the smaller the stable region. All the identified control gains are within the stable region and are distributed close to the node-spiral separation line. Parameter points (p * , d * ) (fastest decay) and (p,d ) (closest point on node-spiral separation line) are indicated by black × and green × markers for each subject. The basic statistical results of the fitted control parameters τ, p and d are collected in Table 1. The best fitting time delays were found to be in the range of 100 ∼ 200 ms, which is in agreement with different estimates in the literature [48][49][50][51] . The identified control gain values also resemble those in the literature for delayed PD feedback models of human quiet standing 50,51 .

Settling time versus robustness to parameter changes
Experiments showed that the fitted control gains are slightly larger than the gains (p * , d * ) corresponding to the fastest decay. An explanation for this is that the sensitivity of the exponential decay rate γ 1 to parameter changes is different at different    sections of the node-spiral separation line. Note that the response is bounded by |θ (t)| θ 0 e γ 1 t , therefore changes in γ 1 are amplified through an exponential function. The change of γ 1 along the node-spiral separation line is shown in Figure 4 for the parameters of Subject 9. Panel a shows the stability diagram with some sample values of γ 1 . It can be seen that γ 1 changes faster in the lower (black) section of the node-spiral separation line than in the upper (green) one. Panels b and c show the change of γ 1 as function of p and d, respectively. The system with control gains selected from the lower branch is more sensitive to changes in p and d than the system corresponding to the upper (green) branch. Actually, the parameter point (p * , d * ) is infinitely sensitive to negative changes in p and d, which can be seen from the vertical slope of its tangent in panels b and c. When applying ±10% perturbation in the gains (p * , d * ), then, in the worst case, γ 1 changes from γ * = −2.98 to −1.84. The same ±10% perturbation on the gains (p,d), in the worst case, results in changes of γ 1 from −2.59 to −2.38. Hence, when perturbations/uncertainties are present in the sensory perception, then it is more beneficial to select the control gains form the upper (green) branch of the node-spiral separation line than from the lower (black) one.
Based on this observation, three measures are introduced to describe the response in terms of the control gains. Let the mean of the fitted control gains be denoted by (p,d). First, the oscillatory nature of the response related to hypothesis H1 can be described as the distance between the point (p,d) and closest point (p,d) of the node-spiral separation line (see Figure 1). Point (p,d) is defined such that (p −p) 2 + (d −d) 2 has to be minimal (for this, the control gains have to be normalized by introducing the dimensionless timet = t/τ). Second, the decay of the response (i.e., settling time) related to hypothesis H2 can be characterized by the distance between the fitted parameter point (p,d) and the parameter point (p * , d * ) associated with the fastest decay. Third, the combined concept of fast decay and robustness to uncertainties in the control gains related to hypothesis H3 can be characterized by the relationsp > p * andd > d * .
The above discussion implies that the goal function might not be to achieve (p, d) = (p * , d * ) but to tune the control gains to (p, d) = (p,d), which guaranties a lower limit to γ 1 even in the case of perturbations of the control gains. This concept minimizes the settling time (e.g., maximizes the magnitude of γ 1 ) while at the same time preserves robustness to static perturbations in the control gains (p, d).

5/14
Hypothesis 1: node-spiral separation line -accepted Control gains associated with oscillatory and non-oscillatory responses are separated by the node-spiral separation line. Hypothesis 1 can be tested using the distance between the identified mean control gains (p,d) and the closest point (p,d) on the node-spiral separation line. The corresponding one sample t-test (ttest(p,p) and ttest(d,d)) are shown in Table 2. Results show that Hypothesis H1 related to the location of the proportional gain p is accepted for most of the subjects (7 out of the 10) and also for the overall data. H1 related to the location of the differential gains d is accepted for all the subjects. This observation confirms that control gains are tuned to be close to the node-spiral separation line.

Hypothesis 2: fastest decay of the response -rejected
Statistical analysis was performed in order to check whether the fitted control gain pairs (p,d) correspond to the point (p * , d * ) associated with the fastest decay. The results of one sample t-test (ttest(p, p * ) and ttest(d, d * )) are shown in Table 2. In case of most of the subjects and also in case of the overall data, Hypothesis 2 was rejected. The fitted control gain pairs are not distributed around the control gain pairs that yields the fastest decay. An explanation to this observation is that the parameter point (p * , d * ) is infinitely sensitive to small changes in the control gains, hence to the small changes in the perceived sensory feedback.

Hypothesis 3: fast decay of the response with robustness -accepted
The fitted control gains was shown to be close to the point (p,d) on the node-spiral separation line in H1. Now it is to be checked whether (p,d) lies on the upper or on the lower branch of the separation line. For 9 out of the 10 subjects, it was observed thatp > p * andd > d * , hence (p,d) lies on the upper (green) branch. This confirms that the control gains are indeed tuned to the more robust section of the node-spiral separation line, hence H3 is accepted. Therefore, it is a plausible assumption that the control gains are tuned to achieve fast decay (γ-stability) but at the same time allow some variations in the gains or in the sensory feedback.

No difference between PA and AP parameters
The difference between the fitted control gains obtained for the PA and the AP trials was analyzed using an unpaired (twosample) t-test assuming unequal variances. No significant difference was found (the P-value for p was v p = 0.941 while for d was v p = 0.167). Hence, both AP and PA balancing process can be modeled by delayed feedback with control gains tuned close to (p,d).

Variation of the fitted gains over the trials -no learning
The variation of the fitted control gains over the 10-trial series is shown in Figure 5 in order to check whether the control gains are coherently tuned towards certain region of the plane (p, d). There is no clear trend in the change of the fitted parameters (either in the mean or in the variations), which suggests that learning process was not present during the trials. Hence, reacting to perturbation during standing still can be considered as an already learnt and acquired feedback mechanism.

Discussion
The results confirm that recovery of quiet standing after a sudden perturbation can well be described by a delayed state feedback mechanism described by three parameters: the reaction delay τ, the proportional and the derivative control gains p and d, respectively (see Methods for validation). While the reaction delay is an inherent feature of the control process, the control gains can be tuned to improve performance. Parameter fitting shows that the control gains are tuned such that the response is non-oscillatory (H1) and result in fast recovery even in case of parameter uncertainties related to the perception of the angular position and the angular velocity (H3). Hence, the typical response is fast and non-oscillatory that resembles the dynamics of a critically damped nondelayed mechanical system 47,52 . This feature of the response is also in agreement with delayed models of human standing in the sense that the control gains are located in the lower left regions of the stability region 43,48,53 .
No significant changes were observed on the parameters over the trials implying that the identified feedback mechanism has been already learned and practiced before during the activities of daily living. The effect of learning process is more pronounced when a new and unknown task is to be performed, e.g, ball-and-beam balancing 45 , balancing on balance board 54 ,

7/14
beam walking 55 or combined quiet standing and stick balancing 28 . A question arises whether practice or other techniques can be developed to further improve the performance against sudden perturbations.
It shall be mentioned that other than delayed PD feedback is also a possible concept to describe the stabilization process. Several ideas are available in the literature, such as intermittent feedback where intermittency can be either a control logic [25][26][27][28]56 or an inherent consequence of the sensory system 20, 30 , acceleration feedback 33,57 or predictor feedback 40,44 , just to mention a few. An advantage of the delayed PD feedback model is that while it is widely used in the literature 27,32,34,48,51 , it accounts with the two most important features of neural motor control: (1) actuation is performed based on perceived sensory signals; and (2) there is a reaction time delay between sensory perception and action.
It is important to note that the tilt angle θ observed for human postural sway and the corresponding dead zones (≈ 0.05 − 0.08 • ) are much smaller than the tilt angles used in our experiments (≈ 3 • ). Although the complexities of the dynamics of postural sway provide insights into the nature of the neural control of balance during quiet standing, it is more likely that the responses to larger perturbations are more relevant for understanding the etiology of falls. Thus, we anticipate that the observation that human balance control responses to large perturbations resembels that of a critically damped, non-delayed mechanical system will have important implications. At the very least, the mathematical analysis is greatly simplified.

Mathematical model
An inverted pendulum model for human standing balance is shown in Figure 1a. Briefly the human body is modeled as a homogeneous rod of mass, m, pivoted on a joint A 31,33 . The equation of motion takes the form where h denotes the distance between the center of gravity and the ankle joint A, g is the acceleration due to gravity, J A is the moment of inertia with respect to point A, and θ is the general coordinate which describes the angular position of the body with respect to vertical. There are two types of torques which interact to stabilize the upright position. First, there is the passive stiffness of the ankle related to the mechanical properties of the foot, Achilles tendon and aponeurosis. The contribution of these forces to balance is modeled by a torsional spring of stiffness k t 20, 22, 31 . The intrinsic mechanical stiffness of the ankle is not sufficient to maintain stability during quiet standing and contractions of parallel connected calf muscles are required. Thus active muscle contractions produce feedback-driven torques, Q(t), which act across the ankle joints. The proportional-derivative feedback control had the form where τ is the reaction time delay, and P and D are, respectively, the proportional and derivative gains. Substitution into (1) and linearization about the θ = 0 upright vertical position gives the delay-differential equation is a system parameter and are normalized control gains. The stability of (3) can be determined using the D-subdvision method 41 . Substitution of the exponential trial solution θ (t) = Be λ t into (3) gives the characteristic function The characteristic equation D(λ ) = 0 has infinitely many complex roots, λ i (i = 1, 2, . . . , ∞), which are called characteristic exponents. The system is stable, i.e, the solution θ (t) converges to 0, if Re(λ i ) < 0 for all i. Taking λ = γ ± iω and setting γ = 0, the equation D(λ ) = 0 gives the transition curves if ω = 0 : p(ω) = (ω 2 + a) cos(ωτ), d(ω) = ω 2 + a ω sin(ωτ)

8/14
These parametric curves delimit the D-shaped stability region shown in Figure 1b. In order to characterize the response associated with different parameter pairs (p, d), the general solution B i e Re(λ i t) cos(Im(λ i t)) + i sin(Im(λ i t)) has to be analysed, where B i is the complex amplitude corresponding to λ i . The values of B i 's are determined by the initial functions (initial perturbations) during t ∈ [−tau, 0], but the stability is independent of these parameters. It is assumed that the characteristic exponents are ordered such that Re(λ 1 ) ≥ Re(λ 2 ) ≥ Re(λ 3 ) ≥ . . . . The dynamics of the response is determined by the dominant (rightmost) characteristic exponent λ 1 . While the real part Re(λ 1 ) corresponds to the decay of the response (settling time), the imaginary part Im(λ 1 ) gives the oscillation frequency. In mathematical terminology, γ 1 = Re(λ 1 ) < 0 is called exponential decay rate and the system is said to be γ-stable if γ 1 ≤ γ < 0 42 . Figure 1 shows the dynamic behaviour of (3). The stable parameter region (where γ 1 < 0) is bounded by black thick curve. Control gains out of the stable region results in either increasing oscillations or exponential growth, hence, in both cases, falling. Within the stable region, thin curves represent different contour lines of γ-stability. The larger the magnitude of γ 1 the shorter the settling time. The fastest response is obtained when (p, d) = (p * , d * ). Thus, p * and d * can be considered as optimal gains with respect to settling time. This concept can be associated to the terminology of critical damping of nondelayed models, which plays an important role in modelling the response to sudden perturbations during standing still 47,48,52 .
The stable region can be separated into two parts based on the imaginary part of the dominant root λ 1 . Darker shaded region to the left from the black-green dashed line is associated with Im(λ 1 ) = 0. In this case, the dominant solution component B 1 e λ 1 t is non-oscillatory. Lighter shaded region to the right of the dashed line is associated with Im(λ 1 ) > 0 and the corresponding solution component reads B 1 e λ 1 t +B 1 eλ 1 t , which is oscillatory with angular frequency Im(λ 1 ) (hereλ 1 andB 1 are the complex conjugate of λ 1 and B 1 , respectively). The black-green dashed line is called node-spiral separation line since it separates node type and spiral type solutions 45 . Note that the point (p * , d * ) corresponding to the fastest response lies on the node-spiral separation line. The node-spiral separation line can be divided into two parts based on the location of the dominant (rightmost) roots. In the lower branch (black dashed line) in Figure 1b, the rightmost characteristic root is real and has a multiplicity of 2.
In the upper branch (green dashed section), a real and a complex pair of characteristic roots coexists with the same real part. At parameter point (p * , d * ), the rightmost root is real (λ 1 = γ * 1 ) and has a multiplicity of 3. Besides the actual values of the exponential decay rate γ 1 , its robustness to changes in the control parameters is also an important feature of the control process. ε p relative error in p and ε d relative error in d alter the control force as hence this perturbation can also be implemented as perturbation in the sensory perception of θ andθ with the same relative error ε p and ε d , while the gains p and d are constant. This suggest a sensitivity analysis of γ 1 to changes in p and d.

Participants
We carried out the experiments with 10 subjects (8 males, 2 females) whose parameters and related statistical data are shown in Table 3. The subjects had no self-reported medical conditions which could affect their ability to perform the required tasks. The research was carried out in accordance with relevant guidelines and regulations following the principles of the Declaration of Helsinki. All subjects provided written informed consent for the procedures, signed a General Data Protection Regulation (GDPR) form and were given the opportunity to withdraw from the study at any time. The research project and the study protocol was approved by the Faculty of Mechanical Engineering, Budapest University of Technology and Economics.

Procedure
The concept of the measurements is shown in Figure 2. We perturbed standing balance by the unexpected release of a resisting force 8 . While standing comfortably the subject resists a horizontal force, F H provided by hanging a weight via a rope that was connected to the subject by a body harness. A constant force F H was applied either in the PA or in the AP direction (Figure 2a and 2b, respectively). Under these conditions the subject's preferred standing position was slightly tilted in order to resist the applied force. The force was released manually using a bolt mechanism at an unexpected moment. This causes an initial sway in the direction opposite to the released force. After some transients, subjects found their new vertical equilibrium (normal posture) and they kept on standing in this position for t s = 15 s. Subjects were instructed to keep their hip and knee joints in a constant extended position and to recover their balance without flexing their knees or hips or moving their arms. The applied force F H shown in Table 3 was the largest one that the subject was able to resist without any difficulties. Each subject performed 10 trials with F H applied in one direction (either PA or AP) followed by 10 trials in which F H was applied in the other direction. In order to prevent carry order effects, the direction of F H for the first 10 trials was chosen randomly. The time instant when the weight was released is t 0 . Maximum excursion was reached at time instant t = t 1 , normal posture was recovered at t = t 2 , after that subjects kept on standing quietly for t s = t 3 − t 2 = 15 s and the balancing trial was terminated after time instant t = t 3 . Response was evaluated and parameter fitting was performed over the period t ∈ [t 1 ,t 1 + t p ] with t p = 10 s. The time signal over the period t ∈ [t 2 ,t 3 ] was used to calibrate the normal posture such that the mean value of the tilt angle θ (t) during this period was zero.

Data collection
A high-speed motion capture system (8 synchronized OptiTrack Prime13 cameras, 120 Hz) was used to measure the three dimensional position (x i (t), y i (t), z i (t)) of spherical reflective markers (diameter 16 mm) where the subscript i refers to the location of the markers: i = 0, 1, 2, 3, 4 respectively stand for the ankle, knee, hip, shoulder, and head. Projecting all movements to the anterior-posterior plane (x, z), the absolute tilt angles can be calculated as where i = 1, 2, 3, 4 respectively indicate ankle-knee, ankle-hip, ankle-shoulder and ankle-head angles.

Validity of the single inverted pendulum model
Tilt angles θ i (i = 1, 2, 3, 4, ankle-knee, ankle-hip, ankle-shoulder and ankle-head angles) were calculated based on the marker positions located on the ankle, knee, hip, shoulder and head, respectively. Measurements show that most of the corrective motion happen at the ankle joint as the participants were instructed to keep the knee and hip joints fixed. The angle across the knee joint was found to be small, i.e., θ 1 − θ 2 ≈ 0, and was further reduced at perturbation onset by contraction of the rectus femoris muscle 46 . Hence, θ 2 ≈ θ 1 provides a good measure of the ankle across the ankle joint. The distribution of ankle-hip angle θ 2 and the difference between the ankle-shoulder and the ankle-hip angles (θ 3 − θ 1 ) are shown in Figure 6. In average, about 75% of the corrective motion takes place at the ankle joint and only ∼ 25% at the hip joint. Correlation between the angle across the ankle joint (θ 2 ) and the ankle-hip angle (θ 3 ) for all subjects were above 0.9, which reflects that corrections at the ankle and hip joints are performed in phase. These observations suggest that a single inverted pendulum model can be used to capture the main characteristics of postural sway 31,34,35,51 . The tilt angle associated with the single pendulum model of the body was calculated as the average of the ankle-hip, ankle-shoulder and ankle-head angles: θ (t) = 1 3 (θ 2 (t) + θ 3 (t) + θ 4 (t)) .