 Research
 Open Access
 Published:
Frequencydependent force direction elucidates neural control of balance
Journal of NeuroEngineering and Rehabilitation volume 18, Article number: 145 (2021)
Abstract
Background
Maintaining upright posture is an unstable task that requires sophisticated neuromuscular control. Humans use foot–ground interaction forces, characterized by point of application, magnitude, and direction to manage body accelerations. When analyzing the directions of the ground reaction forces of standing humans in the frequency domain, previous work found a consistent pattern in different frequency bands. To test whether this frequencydependent behavior provided a distinctive signature of neural control or was a necessary consequence of biomechanics, this study simulated quiet standing and compared the results with human subject data.
Methods
Aiming to develop the simplest competent and neuromechanically justifiable dynamic model that could account for the pattern observed across multiple subjects, we first explored the minimum number of degrees of freedom required for the model. Then, we applied a wellestablished optimal control method that was parameterized to maximize physiologicallyrelevant insight to stabilize the balancing model.
Results
If a standing human was modeled as a single inverted pendulum, no controller could reproduce the experimentally observed pattern. The simplest competent model that approximated a standing human was a double inverted pendulum with torqueactuated ankle and hip joints. A range of controller parameters could stabilize this model and reproduce the general trend observed in experimental data; this result seems to indicate a biomechanical constraint and not a consequence of control. However, details of the frequencydependent pattern varied substantially across tested control parameter values. The set of parameters that best reproduced the human experimental results suggests that the control strategy employed by human subjects to maintain quiet standing was best described by minimal control effort with an emphasis on ankle torque.
Conclusions
The findings suggest that the frequencydependent pattern of ground reaction forces observed in quiet standing conveys quantitative information about human control strategies. This study’s method might be extended to investigate human neural control strategies in different contexts of balance, such as with an assistive device or in neurologically impaired subjects.
Background
Controlling balance during standing and walking is a fundamental necessity for human mobility. Although maintaining upright posture involves little overt movement, its inherently unstable nature poses an interesting sensorimotor control problem [1,2,3,4].
While many recent studies have investigated balance by applying perturbations to the individual [1, 5,6,7], it is also important to understand how humans maintain their balance without external perturbations, i.e., during quiet standing. In particular, the center of pressure and the fluctuations of the center of mass have been commonly used to evaluate balance performance during quiet standing [8, 9]. However, studying the center of mass and/or the center of pressure trajectories alone is insufficient to describe the complex dynamics and control of the multisegmented human body. Insights can be gained by investigating how humans use the direction of their foot–ground interaction force, which is the outcome of a complex sensorimotor control process that involves timed muscle activity, biomechanical constraints, and sensory feedback from multiple pathways. Importantly, the ground reaction forces directly contribute to the centroidal dynamics of the human body [10]. The orientation of the ground reaction force vector and where its lineofaction lies relative to the center of mass may give further insight into how human subjects control the translational motion of the center of mass and net angular motion of the body.
Recently, Gruben and colleagues suggested a new method to study the relation between the orientation of the ground reaction force vector and the center of pressure in human subjects during quiet standing [11, 12]. Net ground reaction forces at different times, which have different orientations and points of application (centers of pressure), intersect at some point in space. The authors defined this point as the intersection point and examined its relation to the center of mass of the standing individual. Because the height of the intersection point relative to the center of mass determines the translational and angular components of centroidal accelerations, it provides a compact geometric representation that is useful for understanding the dynamics and control of human standing balance. When analyzing the force vectors in the frequency domain, this previous study [11] found that the vertical position of the intersection point exhibited a consistent pattern across subjects. With this observation, the authors suggested that the frequencydependent intersection point characterizes the neural controller of human balance. However, the biomechanics of upright posture might account for some of the variation of intersection point height across different frequency bands. Therefore, this study aimed to elucidate the extent to which the frequencydependent variation of the intersection point could be attributed to neural control strategies or to biomechanics.
To this end, the first objective was to develop the simplest competent and neuromechanically justifiable dynamic model that could account for the consistent pattern observed across multiple subjects [11]. Second, we examined the hypothesis that the neural control strategies in standing balance would economize control effort [6]. To test this hypothesis, we took advantage of the linear quadratic regulator (LQR) [13], a wellestablished optimal control method, that enabled a systematic search of physiologicallyplausible controller parameters [2, 14,15,16].
Methods
Human experiment
Experimental procedure
In the previous study [11], ten unimpaired, young participants (24.2 ± 10.3 years) were asked to stand quietly while viewing a mark at head height 1 m away. Each participant completed one 50 s trial standing on a 6axis forceplate measuring at 1000 Hz. The analysis was confined to the sagittal plane. The subjects’ average mass and height were 71 kg and 1.75 m, respectively.
Intersection point
The intersection point was defined as a point in space where the net ground reaction force vectors at adjacent timesteps intersect [11], as illustrated in Fig. 1. The intersection point is a geometric representation of the relation between the ground reaction force and the center of pressure. This point was originally identified with the goal to understand how humans maintain balance during walking [17, 18]; Gruben and colleagues were the first to apply it to understand the mechanics of standing balance [11].
Assuming subtle movements of the body and small variations in ground reaction force orientations, the orientation of the ground reaction force (\(\theta _F\)) can be approximated as
where \(F_x\) and \(F_z\) are the horizontal and vertical components of the net ground reaction force, respectively.
The height of the intersection point (IP) of two forces at adjacent times (\(F(t), F(t+\delta t)\)) is
where CoP(t) is the center of pressure at time t.
Frequencydependence of the intersection point
Investigating the system response in the frequency domain often yields insight into the structure of a dynamic system. To parse the time series into frequency bands, a Hamming window with the length of the entire data set was first applied to both \(\theta _F\) and CoP signals. \(\theta _F\) and CoP signals were then bandpassfiltered (zerolag, 2ndorder Butterworth) and parsed into bands of 0.2 Hz width centered on frequencies from 0.5 to 7.9 Hz (38 nominally nonoverlapping bands). Finally, the principal eigenvector of the bestfit covariance matrix of \(\theta _F\) plotted against CoP (both signals detrended to have zeromean) was extracted for each band. Its slope is equivalent to the inverse of the intersection point, as illustrated in Fig. 2. Assuming small variation between the forces,
the lowerorder component of the intersection point height (IP) can be approximated as
and rearranging (2) results in
Gruben and colleagues [11] found that the vertical position of the intersection point exhibited a consistent pattern across subjects: it was above the center of mass at low frequencies and decreased as frequency increased, reaching an asymptote below the center of mass at higher frequencies as shown in Fig. 4a in the Results section.
Simulation
Single inverted pendulum model
We first investigated whether a single inverted pendulum, which is a widely accepted model for human quiet standing [4], could reproduce the experimental observations. Theoretical analysis showed that the model could not adequately reproduce the experimental observation in [11], because the intersection point height of the single inverted pendulum was always above the center of mass (Appendix 1). Hence, a multidegreeoffreedom model was required.
Double inverted pendulum model
The double inverted pendulum model that was used to simulate a multisegmented human body is illustrated in Fig. 3. The lumped model parameters summarized in Table 1 used the anthropometric distribution of male subjects in the sagittal plane [19] based on the average height and weight of the subjects from [11]. Any mass and length below the ankle was neglected, as the simulation assumed the ankle to be a pin joint. The center of mass positions were measured with respect to the ankle joint for link 1 and the hip joint for link 2. The moments of inertia were calculated about the center of mass of each link.
The equations of motion of the double inverted pendulum were
where \({\mathbf {M}}({\mathbf {q}}) \in {\mathbb {R}}^{2\times 2}\) is the inertia matrix, \({\mathbf {C}}({\mathbf {q}},\dot{\mathbf{q}}) \in {\mathbb {R}}^{2\times 2}\) contains the Coriolis and centrifugal terms, \({\mathbf {G}}({\mathbf {q}}) \in {\mathbb {R}}^{2\times 1}\) are gravitational torques, and \({\varvec{\tau }}=[\tau _1, \tau _2]^T \in {\mathbb {R}}^{2\times 1}\) is the joint torque vector (see Appendix 2 for full symbolic inertia, centrifugal, and gravitational matrices). Generalized coordinates are \({\mathbf {q}} = [q_1, q_2]^T \in {\mathbb {R}}^{2\times 1}\) as defined in Fig. 3. These variables represent the sagittal plane angular displacements of the ankle and hip joints, respectively.
Defining the state vector as \({\mathbf {x}} = [{\mathbf {q}}^T, \dot{\mathbf{q}}^T]^T\), (4) can be rewritten in statedetermined form as
The internal perturbations that cause persistent sway in quiet standing were simulated by additive noise
where \({\varvec{\tau }}_{\mathrm{{ctl}}} = [\tau _{\mathrm{{ctl}},1},\tau _{\mathrm{{ctl}},2}]^T\) are the ankle and hip torques that stabilize the body. In this study, we assumed the noise \({\mathbf {w}} \in {\mathbb {R}}^{2\times 1}\) was white, mutually uncorrelated, and that it followed a zeromean Gaussian distribution with covariance matrix \(E\{\mathbf {ww}^T\}=diag\{\sigma _1^2, \sigma _2^2\}\). The relative strength of the two noise sources was defined as \(\sigma _r =\sigma _1/\sigma _2\), where \(\sigma _1\) and \(\sigma _2\) are the noise at the ankle and hip, respectively.
\(F_x\) and \(F_z\), the horizontal and vertical components of the ground reaction force, were computed as follows
where \(m=m_1+m_2\) is the total mass of the body, \(\ddot{r}_{CoM,x}\) and \(\ddot{r}_{CoM,z}\) are the horizontal and vertical components of the center of mass acceleration. The center of pressure was then computed as
Linear quadratic regulator
This study used a nonlinear model with a linear controller. Hence, the nonlinear equations of motion (5) were first linearized about the upright balancing posture at rest (\({\mathbf {x}}_* = {\mathbf {0}}\) and \({\varvec{\tau }}_* = {\mathbf {0}}\)) as follows
where \(\bar{{\mathbf {x}}} = {\mathbf {x}}  {\mathbf {x}}_*\), \({\varvec{{\bar{\tau }}}} = {\varvec{\tau }}_{\mathrm{{ctl}}}  {\varvec{\tau }}_*\), and \({\mathbf {A}}_{lin}\) and \({\mathbf {B}}_{lin}\) are linearized state and input matrices, respectively (see Appendix 3 for the linearized statespace matrices).
As normal human standing is evidently stable in the upright position, the LQR method was chosen as it guarantees a stable closedloop system^{Footnote 1}. The LQR is an optimal linear statefeedback controller that minimizes the quadratic cost function
to determine control torques
where \({\mathbf {K}}_{LQR}\) is the optimal control gain matrix found via the LQR procedure. The matrices \({\mathbf {Q}}\) and \({\mathbf {R}}\) in (9) weight the state and input deviations from zero.
We parameterized the input weighting matrix as
to facilitate exploration of two important features: the relative cost between state deviation and control effort, determined by the parameter \(\alpha\), and the relative magnitude of hip and ankle effort, determined by the parameter \(\beta\).
When \(\alpha\) is large, control effort is reduced to a minimum value required for stability. Thus, with this choice of \(\alpha\), the need to add joint torque limits to the model was eliminated. Additionally, with large \(\alpha\), the resulting closedloop system has a welldefined behavior (placing its poles at the mirror images of the unstable openloop poles) that is independent of the state weighting matrix \({\mathbf {Q}}\). To evaluate the working hypothesis that humans economize effort, the minimaleffort solution was of interest. Consequently, the choice of the state weighting matrix was not critical, and \({\mathbf {Q}} = I_4\), the identity matrix with dimension 4, was chosen to equally penalize each state’s deviation from equilibrium.
When \(\beta > 1\), the ankle torque is penalized more heavily than the hip, and vice versa when \(\beta < 1\). Since the LQR controller minimizes a quadratic cost function (9) to achieve stability, only the symmetric components of \({\mathbf {R}}\) affect the result. The diagonal values of \({\mathbf {R}}\) were selected such that the size of the matrix (i.e., the product of its eigenvalues) was always equal to 1 and only the components’ ratio affected the results. This choice of parameters also allowed for conclusions to be drawn about the relative penalty on the ankle and hip joints.
Simulation protocol
The simulation was conducted using semiimplicit Euler integration. The initial condition was set to \({\mathbf {x}}_0 = [0, 0, 0, 0]^T\). Replicating the experimental protocol of [11], each simulation was run for 50 s at 1000 Hz. All simulations were conducted in MATLAB 2020a (Mathworks, Natick MA).
To observe the effect of altering the LQR parameters on the frequencydependence of the intersection point and to find the simplest model that could reproduce the human data, various parameters were tested using the following procedure. First, the parameter that weights the relative cost of the control input, \(\alpha\), was set to a large value to ensure minimal control (\(\alpha > 10^4\)). This design choice effectively reduced the number of parameters to two (\(\beta\) and \(\sigma _r\)) as the essential intersection point frequencydependence (above the center of mass at low frequencies, below at high frequencies) varied little when \(\alpha\) was sufficiently large. Then the noise ratio, \(\sigma _r\), was adjusted to produce the best fit at high frequencies while setting \(\beta = 1\). Lastly, \(\beta\) was varied to produce the best fit in the frequency range where the intersection point height was approximately equal to the center of mass height. At the same time, it was ensured that the asymptotic behavior and the fit at high frequencies were maintained.
40 trials were conducted for each tested parameter set to enable statistical analysis of the simulated dependence of the intersection point height on frequency.
Comparison of simulation and human experimental results
When determining the goodness of fit across different model parameter conditions, the average difference of the simulated data and the human subject data from [11] was computed over selected frequency bands by
where \({\text {Human Data}}_i\) is the median of the intersection point height as a fraction of the center of mass height reported by [11] at each frequency band; \(\text {Simulation Data}_i\) is the average intersection point height as a fraction of the center of mass height across 40 trials of the simulation data in a given frequency band; \(N_{band}\) is the number of frequency bands for which the difference in the data was computed. Because balance is characterized by only small motions, a constant center of mass height was assumed.
To identify the onset of the highfrequency asymptote, the human data were fit to an exponential function. The bestfit decay constant was \(T\cong 1\) Hz. Assuming the curve reached its asymptote at frequency \(\cong 3T\), the asymptote started at 3 Hz. Therefore, the difference between the simulated and experimental asymptote was evaluated at frequencies from 3 to 8 Hz (\(N_{band} = 25\)). To evaluate the effect of different controller parameters on the frequency range in which the experimentally observed intersection point height crossed the center of mass height, the average difference between simulation and human data was evaluated at frequencies from 1.2 to 2.6 Hz (\(N_{band} = 7\)). This range encompassed the frequencies in which the observed intersection point height was not statistically different from the center of mass height in [11]. Onesample ttests were used to evaluate the difference between the center of mass height and the simulated mean intersection point height. The 95% confidence interval of the mean of the difference was computed as well.
Results
Minimum required model complexity
Theoretical analysis showed that a singledegreeoffreedom inverted pendulum model could not reproduce the experimental observation in [11]. The intersection point height of a single inverted pendulum model was always above the center of mass for any selection of parameters (Appendix 1). Hence, we proceeded with a double inverted pendulum, i.e. with two degrees of freedom, to approximate the multisegmented human body.
Bestfit model parameter set
The simulated center of mass height did not deviate far from 0.97 m, the height of the center of mass when perfectly upright, justifying the assumption of small angular displacement. In what follows, the center of mass height was therefore assumed to be constant.
The simulated frequencydependent intersection point response for the parameter set, \(\alpha = 10^6\), \(\beta = 0.3\), \(\sigma _r = 0.9\), best matched the human subject data from [11] as shown in Fig. 4. Both simulations and human experimental results show that the intersection point height crossed the center of mass height in similar frequency bands (1.2–2.6 Hz) and had similar asymptotes at higher frequencies. The difference compared to human data for this parameter set was \(0.101 \pm 0.040\) in the 1.2–2.6 Hz range and \(0.011 \pm 0.019\) in the 3–8 Hz range (both within the 95% confidence interval). While the average mass and height in our simulations were taken from [11] to afford best comparison with the human data in Fig. 4c, varying the mass and height parameters over physiologically plausible values did not significantly affect the results.
Varying model parameters
Varying the simulation parameters affected both the frequency at which the intersection point crossed the height of the center of mass and the highfrequency asymptote. The effect of changing parameter values is presented in Fig. 5a–c. The differences between simulation and human data for certain parameter sets are shown in Fig. 5d, e.
Effect of \(\alpha\)
As shown in Fig. 5a, when \(\alpha\), the weighting of control effort relative to state deviation, was increased, the intersection point crossed the center of mass at lower frequencies. For example, when \(\alpha\) was varied from \(10^{4}\) to \(10^6\), the frequency at which the intersection point crossed over the center of mass moved from 3.9 to 1.5 Hz. As expected from theory, when \(\alpha\) was relatively large (\(\alpha > 10^4\)), there was little effect of varying its value on the difference between human and simulation data for different model parameter sets, as shown in Fig. 5d, e.
Effect of \(\beta\)
As shown in Fig. 5b, when \(\beta\) was decreased, i.e. when hip control was penalized more than ankle control, the intersection point crossed the center of mass at higher frequencies. For example, when \(\beta\) was varied from 1 to 0.3, the frequency at which the intersection point crossed over the center of mass moved from 1.1 to 1.5 Hz. In Fig. 5d, \(\beta = 0.2\) was shown to be the parameter with the smallest difference (0.024) in the 1.2–2.6 Hz range when \(\alpha = 10^6\). However, both the selection of \(\beta = 0.2\) and \(\beta = 0.1\) sacrificed the highfrequency fit, increasing the absolute value of the difference in that range by 0.102 and 0.292, respectively, compared to \(\beta = 0.3\) when \(\alpha = 10^6\). As \(\beta\) deviated from \(\beta = 0.3\), the absolute value of the difference in the 1.2–2.6 Hz range increased by 0.181 when \(\beta = 1\) and \(\alpha = 10^6\).
Effect of \(\sigma _r\)
Adjusting \(\sigma _r\) shifted the high frequency asymptote (3–8 Hz) of the intersection point, as shown in Fig. 5c. Variation of the highfrequency asymptote of the intersection point height was predicted by the analysis presented in Appendix 4. Here, the two extremes, zero noise in the ankle (\(\sigma _r=0\)) and the hip (\(\sigma _r=\infty\)), provided lower and upper bounds for the high frequency asymptote. When compared to the bestfit height of the intersection point at high frequencies, the asymptote was 55% higher when \(\sigma _r = 2\) (more noise in the ankle) and 30% lower when \(\sigma _r = 0.5\) (more noise in the hip). In Fig. 5e, \(\sigma _r = 0.9\) is shown to be the bestfit parameter with the smallest difference value, at −0.011 when \(\alpha = 10^6\). As \(\sigma _r\) deviated from the bestfit value, the difference increased to 0.245 when \(\sigma _r = 2\) and to −0.154 when \(\sigma _r = 0.5\), when \(\alpha = 10^6\).
Discussion
This study analyzed a deliberately simplified model of human quiet standing with a stabilizing linear optimal controller to better understand the origin of the frequencydependent intersection point reported by Gruben and colleagues [11].
The simplest competent model required two degrees of freedom (ankle and hip) with a stabilizing controller that used minimal control effort and more ankle torque than hip torque. We successfully identified a narrow range of parameters that provided not only a quantitative reproduction of experimental observations, but also qualitative insight.
Neural control or biomechanics?
Despite the biomechanical constraints that limit the admissible center of mass accelerations and the centers of pressure [14, 20], the ground reaction forces that comply with these constraints are infinite [17, 21]. Beyond the obvious fact that the musculoskeletal system is inherently unstable without a neural controller, we should not expect mechanics alone to determine the intersection point’s frequency dependence. When Gruben and colleagues [11] analyzed the frequency dependence of the intersection point, they observed a consistent trend across multiple subjects and suggested that this consistency was the signature of a neural controller employed by humans during balance.
The results of our simulations replicated the frequency dependence of the intersection point reported for human standing in the sagittal plane—the intersection point was above the center of mass at low frequencies and below the center of mass at high frequencies, as shown in Fig. 4. To understand the general frequencydependent trend of the intersection point, first consider an extreme case at very low frequencies where the twodegreeoffreedom pendulum behaves similar to a single rigid body: its intersection point would be above the center of mass, like that of the single inverted pendulum (Appendix 1). A double inverted pendulum can also be stabilized solely by hip torque, i.e. zero ankle torque. In the latter case, the system would exhibit higher frequency behavior while maintaining the intersection point height to be zero (from (7) and (2) with \(\tau _1=0\)). This indicates that the general trend for the intersection point to be above the center of mass at low frequencies and below at high frequencies may be a consequence of biomechanics, i.e. a double inverted pendulum stabilized about upright posture.
However, biomechanics cannot account for the specific details of the frequency variation of the intersection point height. Somewhere between the lowfrequency and highfrequency regimes, the intersection point must cross from above to below the center of mass height; this particular crossing point is not specified by biomechanics. Similarly, biomechanics does not dictate the asymptote to which the intersection point height converges at high frequencies. In fact, both the intersection point’s asymptote and the frequency at which the intersection point height crossed that of the center of mass varied substantially across the tested parameter values. It was only a small set of model parameters that could replicate human behavior. Therefore, we conclude that the details of the profile of intersection point height with frequency reflect a neural control strategy used by humans during quiet stance.
Physiologicallyplausible bestfit parameters
The main contribution of this work is to deploy a deliberatelysimplified mathematical analysis to elucidate how experimental observations of the frequencydependence of the intersection point may inform the neural control of balance. To conduct this quantitative analysis, the model parameters were systematically varied such that the simulated response of the intersection point's frequencydependence closely replicated human data. To facilitate analysis, we took advantage of the LQR procedure and its wellknown properties.
Selecting \(\alpha = 10^6\) yielded the bestfit result compared to human data, suggesting that a double inverted pendulum model with minimal control effort can successfully account for the frequencydependent intersection point observed in humans. Though it is widely assumed that humans generally minimize effort, supporting evidence during quiet standing has been sparse. Our results show that the observed variation of intersection point height with frequency implies that humans minimize control effort rather than reduce state deviation during quiet standing. This is consistent with the conclusion of a previous study reporting that the nervous system does not exert more control effort than necessary to stabilize upright balance [6].
Long transmission delays in the neural system pose a risk to stability of the balance controller. To account for this, the continuous feedback loop gain must be effectively zero at high frequencies regardless of variations in other model parameters. However, muscle mechanical impedance is not limited in this way; it can respond essentially instantaneously. Behavior in the high frequency range is therefore not likely to depend on neural feedback (defined by \(\alpha\) and \(\beta\)), but instead on neuromuscular impedance and noise (defined by \(\sigma _r\)). Hence, the noise ratio, \(\sigma _r\), was adjusted to fit the highfrequency range before fitting the lowfrequency range with \(\beta\).
Altering the relative noise magnitude in the ankle and the hip torques (\(\sigma _r\)) shifted the highfrequency asymptote of the intersection point height. The simulation result most similar to human experimental data had a 0.9:1 ankletohip noise ratio.
The \(\beta\) value that best described human data in [11] was 0.3 that penalized hip control effort more than ankle control effort. That is, the system is more likely to use the ankle to maintain upright posture than the hip. This is consistent with previous findings that humans primarily use the “ankle strategy” during quiet standing [14, 22,23,24].
Single vs. multijoint model
The observed trend that the intersection point varied with frequency in humans requires multisegment mechanics (Appendix 1). Although the single inverted pendulum model has been widely used to model quiet human standing [1, 4, 26,27,28,29], our finding that a singlesegment model cannot adequately describe quiet standing is also consistent with recent literature [6, 30,31,32,33].
Why no more than two degrees of freedom? It is patently obvious that the standing human body has many more degrees of freedom. However, although adding a knee joint [34] or multiple segments of the spine might more accurately replicate human biomechanics, it is unclear whether this would improve the insight to be gleaned from experimental observations. In fact, as shown in Appendix 4, the twosegment model yields a highfrequency asymptote for the intersection point height that must lie between zero (corresponding to zero noise at the ankle) and below the center of mass height (corresponding to zero noise at the hip). These two extremes bracket the experimental observations reported in [11]. Thus, the twosegment model used in this study was the simplest that could competently reproduce the experimental results observed in [11].
Intersection point: a target variable of control or an emergent consequence?
In this study, the feedback signal was the state error (joint angles and velocities) rather than the intersection point height. Even so, the control model was able to replicate the frequency dependence of the intersection point found in humans. Hence, it appears that the intersection point may be an emergent consequence of stabilization rather than a variable explicitly regulated by the controller. Consistent with this hypothesis, a previous study suggested that the force direction pattern observed in human walking might be an emergent property, rather than a target variable of control [35]. However, further experimentation would be required to test this hypothesis.
Limitations
The simulations conducted in this study assumed simple mechanics. The joint torques in the model are net joint torques that summarize the contributions of various elements, from passive muscle properties to complex neural control. Known features of neuromuscular physiology such as muscle mechanical impedance, neural transmission delay, or sensory noise were omitted. While these features are unquestionably present, our goal was to identify the simplest model competent to reproduce experimental observations. Despite the lack of muscle and nervelevel detail, our simulations were able to articulate subtle differences between control parameters that influence the frequencydependence of the intersection point. Nevertheless, including those neurophysiological features might yield further insight; that is deferred to future work.
This study employed a linear fullstate feedback controller with a constant gain matrix (proportional feedback of angle and angular velocity) even though the central nervous system comprises many nonlinear neural elements. This decision was motivated by the observation that the body generates only small motions about the upright posture, justifying the use of a linearized model to obtain feedback controller gains. This observation also justified the choice of additive noise, as higherorder terms that characterize nonlinear noise processes are negligible. We therefore modeled the noise as white. However, some studies have indicated that biological noise may be better described by ‘pink’ or Brownian noise [25]. Since the lowpass filter property of inertial mechanics suppresses highfrequency components of the spectrum, this noise model proved to be a convenient and viable option.
Finally, the model employed in this study assumed a perfect state estimator. Future studies might assess the effect of including sensory information into the motor controller by employing other control architectures, for instance, using an adaptive [5] or optimal state estimator [14] instead of perfect fullstate measurements.
Another important point to highlight is that we do not presume that the central nervous system actually implements the linear regulator used in our model. The LQR design procedure was simply a tool to generate stabilizing controllers while simultaneously analyzing the influence of factors like the cost of control on balance performance.
Conclusion
This study showed that a double inverted pendulum stabilized by a linear minimaleffort controller could account for the ground reaction force pattern observed in human quiet standing. Numerical simulations also informed the contribution of neural control and biomechanics in generating the pattern observed in human data, i.e. the frequencydependence of the intersection point. The results suggest that the intersection point conveys quantitative information about human balance control strategies.
This study introduced a method to select optimal control and noise parameters that best reproduced human data. This method might be extended to study human neural control strategies in different contexts, e.g., balance in the frontal plane, balance on a beam, balance with and without assistive devices, or in other populations such as aged or neurologically impaired subjects.
Availability of data and materials
Not applicable.
Notes
 1.
To ensure stability, the statespace matrices \({\mathbf {A}}_{lin}\) and \({\mathbf {B}}_{lin}\) must be a controllable pair, the \({\mathbf {Q}}\) matrix must be symmetric positive semidefinite, and the \({\mathbf {R}}\) matrix must be symmetric positive definite.
Abbreviations
 LQR:

Linear quadratic regulator
 CoP:

Center of pressure
 IP:

Intersection point
 CoM:

Center of mass
References
 1.
Peterka RJ. Sensorimotor integration in human postural control. J Neurophysiol. 2002;88(3):1097–118.
 2.
Horak FB. Postural orientation and equilibrium: what do we need to know about neural control of balance to prevent falls? Age Ageing. 2006;35 Suppl 2:ii7–ii11.
 3.
Allen JL, Ting LH. Why is neuromechanical modeling of balance and locomotion so hard? In: Prilutsky BI, Edwards DH, editors. Neuromechanical modeling of posture and locomotion. New York: Springer; 2016. p. 197–223.
 4.
Olsson F, Halvorsen K, Aberg AC. Neuromuscular controller models for quantifying standing balance in older people: A systematic review. IEEE Rev Biomed Eng. 2021.
 5.
van der Kooij H, Jacobs R, Koopman B, van der Helm F. An adaptive model of sensory integration in a dynamic environment applied to human stance control. Biol Cybern. 2001;84(2):103–15.
 6.
Kiemel T, Zhang Y, Jeka JJ. Identification of neural feedback for upright stance in humans: stabilization rather than sway minimization. J Neurosci. 2011;31(42):15144–53.
 7.
Goodworth AD, Peterka RJ. Identifying mechanisms of stance control: a single stimulus multiple output modelfit approach. J Neurosci Methods. 2018;296:44–56.
 8.
Collins JJ, De Luca CJ. Openloop and closedloop control of posture: a randomwalk analysis of centerofpressure trajectories. Exp Brain Res. 1993;95(2):308–18
 9.
Moon J, Pathak P, Kim S, Roh S, Roh C, Shim Y, Ahn J. Shoes with active insoles mitigate declines in balance after fatigue. Sci Rep. 2020;10:1951.
 10.
Orin DE, Goswami A, Lee SH. Centroidal dynamics of a humanoid robot. Auton Robot. 2013;35(2):161–76.
 11.
Boehm WL, Nichols KM, Gruben KG. Frequencydependent contributions of sagittalplane foot force to upright human standing. J Biomech. 2019;83:305–9.
 12.
Yamagata M, Gruben K, Falaki A, Ochs WL, Latash ML. Biomechanics of vertical posture and control with referent joint configurations. J Motor Behav. 2020;53:72–82.
 13.
Sontag ED. Mathematical control theory. New York: Springer; 1998.
 14.
Kuo AD. An optimal control model for analyzing human postural balance. IEEE Trans Biomed Eng. 1995;42(1):87–101.
 15.
Todorov E, Jordan MI. Optimal feedback control as a theory of motor coordination. Nat Neurosci. 2002;5:1226–35.
 16.
Todorov E. Optimality principles in sensorimotor control. Nat Neurosci. 2004;7:907–15.
 17.
Gruben KG, Boehm WL. Force direction pattern stabilizes sagittal plane mechanics of human walking. Hum Mov Sci. 2012;31(3):649–59.
 18.
Maus HM, Lipfert SW, Gross M, Rummel J, Seyfarth A. Upright human gait did not provide a major mechanical challenge for our ancestors. Nat Commun. 2010;1:70.
 19.
de Leva P. Adjustments to ZatsiorskySeluyanov’s segment inertia parameters. J Biomech. 1996;29(9):1223–30.
 20.
Kuo AD, Zajac FE. Human standing posture: multijoint movement strategies based on biomechanical constraints. Prog Brain Res. 1993;97:349–58.
 21.
Duarte M, Sternad D. Complexity of human postural control in young and older adults during prolonged standing. Exp Brain Res. 2008;191(3):265–76.
 22.
Horak FB, Nashner LM. Central programming of postural movements: adaptation to altered supportsurface configurations. J Neurophysiol. 1986;55(6):1369–81.
 23.
Nashner LM, McCollum G. The organization of human postural movements: a formal basis and experimental synthesis. Behav Brain Sci. 1985;8(1):135–50.
 24.
Runge CF, Shupert CL, Horak FB, Zajac FE. Ankle and hip postural strategies defined by joint torques. Gait Posture. 1999;10(2):161–70.
 25.
van der Kooij H, Peterka RJ. Nonlinear stimulusresponse behavior of the human stance control system is predicted by optimization of a system with sensory and motor noise. J Comput Neurosci. 2011;30(3):759–78.
 26.
Morasso P, Cherif A, Zenzeri J. Quiet standing: the Single Inverted Pendulum model is not so bad after all. PLoS One. 2019;14(3):e0213870.
 27.
Morasso PG, Schieppati M. Can muscle stiffness alone stabilize upright standing? J Neurophysiol. 1999;82(3):1622–6.
 28.
Winter DA, Patla AE, Rietdyk S, Ishac MG. Ankle muscle stiffness in the control of balance during quiet standing. J Neurophysiol. 2001;85(6):2630–3.
 29.
Kiemel T, Oie KS, Jeka JJ. Multisensory fusion and the stochastic structure of postural sway. Biol Cybern. 2002;87(4):262–77.
 30.
Rozendaal LA, van Soest AJ. Stabilization of a multisegment model of bipedal standing by local joint control overestimates the required ankle stiffness. Gait Posture. 2008;28(3):525–7.
 31.
Günther M, Grimmer S, Siebert T, Blickhan R. All leg joints contribute to quiet human stance: a mechanical analysis. J Biomech. 2009;42(16):2739–46.
 32.
Pinter IJ, van Swigchem R, van Soest AJ, Rozendaal LA. The dynamics of postural sway cannot be captured using a onesegment inverted pendulum model: a PCA on segment rotations during unperturbed stance. J Neurophysiol. 2008;100(6):3197–208.
 33.
Sasagawa S, Shinya M, Nakazawa K. Interjoint dynamic interaction during constrained human quiet standing examined by induced acceleration analysis. J Neurophysiol. 2014;111(2):313–22.
 34.
Yamamoto A, Sasagawa S, Oba N, Nakazawa K. Behavioral effect of knee joint motion on body’s center of mass during human quiet standing. Gait Posture. 2015;41(1):291–4.
 35.
Müller R, Rode C, Aminiaghdam S, Vielemeyer J, Blickhan R. Force direction patterns promote whole body stability even in hipflexed walking, but not upper body stability in human upright walking. Proc R Soc A Math Phys Eng Sci. 2017;473(2207):20170404.
Acknowledgements
We would like to thank Professor Kreg Gruben for helpful discussion.
Funding
This work was supported in part by NSFCRCNS 1724135, awarded to Neville Hogan and NSFCRCNS1723998, awarded to Dagmar Sternad. Jongwoo Lee was supported by a Samsung Scholarship.
Author information
Affiliations
Contributions
KS performed simulation and data analysis and drafted the manuscript. JL developed simulations and performed mathematical analysis of the intersection point. MR contributed to the frequency analysis of the intersection point. JL, MR, DS, and NH contributed to data and statistical analyses. KS and JL generated figures. All authors edited the final manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix 1: Intersection point of the single inverted pendulum
The intersection point below the center of mass at high frequencies observed in human data cannot be reproduced by a single inverted pendulum model. Consider a single inverted pendulum with mass m, center of mass position from the pivot \(l_c\), moment of inertia about pivot \(j'\), gravitational acceleration g, and actuated by ankle torque \(\tau\). The equation of motion is
where q is the angular displacement of the ankle joint with respect to the upright equilibrium posture, \(\tau _{\mathrm{{ctl}}}\) is the control torque, and w is the additive actuation noise. For small motions typical of quiet standing, linearization of (12) is well justified:
As introduced in (3), the intersection point is defined in terms of the orientation of the force, \(\theta _F\), and the center of pressure, CoP:
Taking the Laplace transform:
where s is a complex variable, Q(s), \(\Theta _F(s)\), and COP(s) are the Laplace transforms of q, \(\theta _F\), and CoP, respectively. Denote \(H(s)=Q(s)/W(s)\), the transfer function from input noise to output motion, where W(s) is the Laplace transform of w. To investigate the intersection point at high frequency, consider \(s=i\Omega\) where \(i^2=1\) and \(\Omega \rightarrow \infty\). As the first term of COP(s) dominates, \(\Theta _F(s)\) and COP(s) have the same phase. Then, the variation of two output variables will be linear at high frequencies and the intersection point can be determined from the ratio of magnitudes of the two outputs:
As \(\Omega \rightarrow \infty\),
Note that the centroidal moment of inertia is \(j=j'mc^2\). Therefore, the intersection point height must be greater than the center of mass height. The single inverted pendulum model cannot explain the intersection point behavior observed in human studies.
Appendix 2: Nonlinear model equations
The equations of motion of the double inverted pendulum are expressed in (4). Note that the choice of generalized coordinates, \(q_1\), \(q_2\), are consistent with the generalized forces (torques) that are applied. The following details each component of the matrices in terms of the variables summarized in Table 1. \(j'_1\) and \(j'_2\) denote moments of inertia taken about ankle and hip joints, respectively. \(l_{c1}\) is the distance from the ankle joint to the center of mass of link 1, and \(l_{c2}\) is the distance from the hip joint to the center of mass of link 2. \(\cos {(q_i)}\) and \(\sin {(q_i)}\) are replaced with \(c_i\) and \(s_i\), respectively. Setting \(\theta _2 = q_1 + q_2\), \(\cos {(\theta _2)}\) and \(\sin {(\theta _2)}\) are replaced with \(c_{\theta _2}\) and \(s_{\theta _2}\), respectively.
The ground reaction forces are computed from the motion of the center of mass, \({\mathbf {r}}_{CoM}({\mathbf {q}})\in {\mathbb {R}}^{2\times 1}\). The acceleration of the center of mass was computed as
\({\mathbf {J}}_{CoM}\in {\mathbb {R}}^{2\times 2}\) is the Jacobian of the center of mass with respect to the joint angles \({\mathbf {q}}\) and \(\mathbf{x} = [{\mathbf{q}},\dot{\mathbf{q}}]^{T}\) is the state vector.
The Jacobian matrix and its derivative are given as follows:
and
where
and \(M_1 = \frac{m_1}{m_1+m_2}\) and \(M_2 = \frac{m_2}{m_1+m_2}\).
Appendix 3: Linearized statespace matrices
Linearizing the equations of motion about the stable upright position, we are left with (8). The statespace matrices are
where
and \({\mathbf {B}} = {\mathbf {I}}_2\).
Appendix 4: Intersection point of the linearized double inverted pendulum
Noting (1) and (7), consider two outputs \({\mathbf {y}}=[y_1, y_2]^T\):
From (13),
Then, linearized output equations can be obtained as
where \({\mathbf {C}}_1 = {\mathbf {J}}_{y_1} {\mathbf {A}}_{lin} {\text{and}} {\mathbf {D}}_1 = {\mathbf {J}}_{y_1}{\mathbf {B}}_{lin}\), evaluated at \(({{\mathbf {x}}, \varvec{\tau })=({\mathbf {x}}^*, \varvec{\tau }^*})\), and \({\mathbf {C}}_2 = {\mathbf {0}} {\text{and}} {\mathbf {D}}_2 = [1, 0]\). With controller \({\varvec{\tau }}=\mathbf {Kx+w}\) as in (6) and (10), the closedloop linear system can be constructed as
The multiinput, multioutput (MIMO) transfer function can be obtained
where \({\mathbf {Y}}(s)\) and \({\mathbf {W}}(s)\) are the Laplace transforms of \({\mathbf {y}}\) and \({\mathbf {w}}\), respectively. The intersection point at each frequency can be obtained by the procedure outlined in the Methods section, as the inverse of the slope of the principal eigenvector. If \(y_1(t)\) and \(y_2(t)\) are harmonic, this procedure is equivalent to finding the slope of the major axis of the ellipsoid that the two signals form.
Assuming two harmonic signals \(y_i(t)\) with magnitude \(\nu _i\) and phase \(\phi _i\) at frequency \(\Omega\),
an implicit formula for the ellipsoid can be written in quadratic form,
where \(\phi =\phi _1\phi _2\). The eigenvector corresponding to the smaller eigenvalue is the major axis and its slope is the inverse of the intersection point as in Fig. 2.
Consider two extreme cases where ankle noise is zero (\(w_1=0\) and \(\sigma _r = 0\)) and hip noise is zero (\(w_2=0\) and \(\sigma _r = \infty\)). For example, when hip noise is zero, substituting \(s = i\Omega\),
and
The intersection point height can be calculated using this method at different frequencies as shown in Fig. 6. Since we are examining a linearized model, nonzero ankle and hip noise responses would be some combination of these two extreme responses.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Shiozawa, K., Lee, J., Russo, M. et al. Frequencydependent force direction elucidates neural control of balance. J NeuroEngineering Rehabil 18, 145 (2021). https://0doiorg.brum.beds.ac.uk/10.1186/s12984021009072
Received:
Accepted:
Published:
DOI: https://0doiorg.brum.beds.ac.uk/10.1186/s12984021009072
Keywords
 Posture and balance
 Inverted pendulum models
 Ground reaction forces
 Neural control
 Biomechanics