Direct estimation of the parameters of a delayed, intermittent activation feedback model of postural sway during quiet standing

Summary

Human postural sway during quiet standing has been characterized as a proportional-integral-derivative controller with intermittent activation. In the model, patterns of sway result from both instantaneous, passive, mechanical resistance and delayed, intermittent resistance signaled by the central nervous system. A Kalman-Filter framework was designed to directly estimate from experimental data the parameters of the model’s stochastic delay differential equations with discrete dynamic switching conditions. Simulations showed that all parameters could be estimated over a variety of possible data-generating configurations with varying degrees of bias and variance depending on their empirical identification. Applications to experimental data reveal distributions of each parameter that correspond well to previous findings, suggesting that many useful, physiological measures may be extracted from sway data. Individuals varied in degree and type of deviation from theoretical expectations, ranging from harmonic oscillation to non-equilibrium Langevin dynamics.

Keywords:

Biology and life sciences – Physical sciences – Engineering and technology – Research and analysis methods – Neuroscience – Cognitive science – Computer and information sciences – Mathematics – Simulation and modeling – Anatomy – Medicine and health sciences – Cognitive neuroscience – Applied mathematics – Algorithms – musculoskeletal system – Body limbs – Signal processing – Systems science – Motor reactions – Postural control – Legs – Control engineering – Ankles – Kalman filter – Control theory – Noise reduction – Differential equations – Numerical analysis – Interpolation

Introduction

Several previous studies have analyzed bodily sway patterns in quiet standing, and a variety of models have been proposed. In this study, we designed and tested a method of directly estimating the parameters of the Asai et al. [1] intermittent feedback control model of posture from experimental data. We begin with a brief review of prior models and the rationale for choosing a model of intermittent postural control (IPC). In the second section, we describe the current model in more detail and explain our framework for the estimation of its parameters. The third section describes simulation studies that tested the estimation capabilities of our framework when the data-generating parameters are known and the model is specified accurately. In the fourth section, we applied the model to experimental data and estimated sampling distributions for each parameter.

Observed trajectories of postural sway have largely been studied as a problem of stochastic behavior, though some studies have focused on its chaotic properties [2]. In this study, we too regarded postural sway as a random process subject to statistical analysis. The center of pressure (COP) on a force plate during quiet standing has been shown to exhibit the features of a bounded, random walk, or correlated noise [3]. Center of mass (COM) is one of the most common metrics of body sway but has to be inferred from other position and force metrics such as the COP [4]. For the small radius in which postural sway occurs, body tilt angle is nearly equivalent to COM and has likewise been used to develop models of posture [1].

Many authors have observed that sway follows a two-frequency oscillation scheme, with fast oscillations of the COP around a drifting center point [3, 5–7]. Collins and De Luca [3] regarded these patterns as a combination of short-term, open-loop system with long-term closed-loop control. Alternatively, the “rambling and trembling” hypothesis suggests that short-term tremors result from corrective, closed-loop feedback that is activated with deviation of the COP from the ground projection of the COM, which is itself allowed to drift [7].

Broadly, more recent debate over the control scheme of human balance has focused on two kinds of models: continuous and intermittent feedback controllers. Continuous control is exerted through a proportional-integral-derivative (PID) or closed-loop system often characterized by a second order linear differential equation, sometimes including delayed proportional and derivative feedback. For instance, Maurer and Peterka [8] tested a PID inverted pendulum model that distinguishes passive, instantaneous feedback from sources such as ankle joint stiffness, from delayed, active feedback from the central nervous system and subsequent muscular response. Others have argued that human postural movement is better described by intermittent feedback mechanisms due to a smaller reliance on process noise, reproduction of cyclical behavior over multiple timescales, more efficient energy expenditure, and greater robustness to disturbances and instability caused by delays in neural signaling [9]. Simulations [10] and reinforcement learning [11] have been used to show that an upright pendulum, taken as a simple model of the standing human body, can exhibit stability and the observed slow oscillation patterns as a result of learned, time-delayed, intermittent feedback.

Intermittent activation models have taken multiple forms. Gawthrop and Wang [12] initially proposed clock-driven muscular feedback, but later considered event-driven models [13]. Event-driven models are generally defined by a combination of stable and unstable manifolds in the phase space of the body’s position or angle. Gawthrop et al. [13] and Eurich and Milton [14] describe the behavior of systems with position-based thresholds that result in two stable equilibria. A model by Bottaro et al. [15] proposes boundary functions of both position and velocity that jointly determine probabilistic bursts of negative feedback. Asai et al. [1] reproduced a commonly observed double power-law structure in the PSD of sway [6] using similar control manifolds but with deterministic rules for sustained feedback activation. Their model requires only a simpler, Gaussian distribution of process noise with a smaller variance as compared to continuous PID models. Nomura et al. [16] showed that the same intermittent activation feedback model is capable of reproducing both chaotic and stochastic patterns that resemble human postural sway as a function of small hemodynamic perturbations, while continuous feedback models cannot.

A common method of estimating the parameters of each model is to simulate data that optimally resemble the experimental data. This is accomplished by varying parameters over iterations of simulation until resulting disparities on a set of key summary statistics have been minimized. Bottaro et al. [15] used the the Root Mean Square (RMS) of both the COP and COM series and each of its derivatives, unimodality of the series histogram, the length of largest oscillations calculated from zero crossings, and the PSD of the COP series. Maurer and Peterka [8] estimated parameters from observed data in a similar manner using mean velocity, RMS distance and velocity, spectral properties such as mean frequency, frequency dispersion, and total power. Asai et al. [1] used the double power law structure of the frequency spectrum as a criterion for the success of their model but did not demonstrate a direct empirical application. To obtain statistical information about estimated parameters, summary statistic methods have been combined with approximate Bayesian inference [17]. This method was used to acquire empirical posterior distributions of five out of the eight parameters of interest [18].

While the simulation approach is flexible for a wide range of model specifications and levels of complexity, it risks overlooking attributes of the data that do not have specific effects on the chosen summary statistics. Bottaro et al. [15] notes, for instance, that “The intermittent nature of the control process cannot be detected by global descriptors of the sway patterns, like the PSD of the COP, because they cannot distinguish between asymptotic and bounded stability”. Furthermore, the amplitude spectrum is invariant to reversal of the signal, giving identical results for potentially different mechanisms of variation. This is problematic when the system includes discontinuous dynamics, such as a sharp impulse followed by a more gradual decay. An alternative approach that better accounts for fine-grained sequential dependence is to estimate the structural parameters directly from the data using Kalman filtering or other iterative techniques. No simulation or descriptive statistics are necessarily used, rather the structural parameters are estimated by minimizing an objective function such as the squared prediction error, or by maximizing the likelihood of the data according to an expected noise distribution. The results obtained by this approach can be sensitive to the exact predictive mechanisms specified in the model, and post-hoc analyses of the estimates can be highly informative about the types and degrees of misspecification. Direct estimation (sometimes called exact estimation by comparison [19]) may be particularly useful when the dynamic structure cannot represented by any descriptive statistics with sufficient specificity. The Asai et al. [1] model of posture may present one such case in that it postulates dependence of the spectral power-law property upon nonlinear, physiological mechanisms of feedback control and their properties. Such properties include the delay in neural signaling, the sensitivity of feedback activation, and the strength of passive versus active corrective forces. Furthermore, the process noise distribution of the Asai et al. [1] model is a Gaussian process and thus accords with the statistical assumptions of the Kalman filter. A drawback of direct estimation is that a misspecified model is not guaranteed to result in any interpretable or accurate parameter estimates if the parameters are highly dependent. If the parameter estimates deviate significantly from their theorized values, we may nonetheless analyze the behaviors they imply and draw general inferences.

The aim of this study was to validate and apply a method of directly estimating parameters for event-driven control with specific focus on the popular intermittent control model by Asai et al. [1]. Validation of this analytic strategy will set a foundation for estimating the parameters of alternative models and more comprehensive comparisons. Following the validation study, we estimated empirical values of each parameter from publicly available COM data [20] and compared our results with theoretically expected values from the literature. We included two previously demonstrated covariates in our analysis, visual feedback and age, to attempt to replicate previous findings as further evidence for the validity of the model.

Model

The intermittent postural control (IPC) model by Asai et al. [1] describes a tension between toppling torque due to gravity and a combination of active and passive resistance mechanisms. Passive resistance is proposed to come from leg stiffness and joint friction and is modeled with instantaneous relations between position, velocity, and acceleration. Active feedback control is proposed to arise from motor responses signaled by the central nervous system and is consequently delayed by about 190-210 ms [21].

The model is provided in terms of body tilt angle (θ) as follows:

The active control terms, f_P and f_D, intermittently respond to θ on a time lag of τ ≈ 200 ms according to the conditions:

The first condition represents a threshold dividing the saddle-type attractor of the toppling acceleration into stable and unstable manifolds. The stable manifold briefly occurs when the tilt angle is moving toward zero, while the unstable manifold is characterized by falling away from zero. The angle of the dividing line is given by the slope parameter α. The second condition describes a radius (r) about the origin within which the tilt angle is too small to be detected or too stable for immediate correction (note that r has conventionally been used to denote the delay time interval in the delay differential equation literature. Here we have preferred τ for that purpose.). By converting the switching threshold slope α into the angle a as α = sin ( a π ) cos ( a π ), we change the upper and lower estimation bounds from [−∞, ∞] to [0, 1]. This way, the parameter a represents the percentage of the phase space, not including the insensitivity radius, for which the active control parameters are non-zero.

The estimable parameters of the SDDE are summarized in Table 1. Many of the parameters have previously been estimated in a variety of ways, sometimes with highly varied results. Tietäväinen et al. [18] used the approximate Bayesian inference [17] with data simulation to estimate P, D, a, τ, and σ. Among these, the method failed to obtain precise distributions for D in both simulations and empirical application. It is also not clear whether fixing the other parameters to uncertain theoretical priors (K = .8, B = 4, and r = .004) results in biased estimates. Direct physiological measurements found the relative resistance to toppling torque at the ankle, K, to be as high as 91% on average [22] when the average magnitude of disturbance is small. Another study estimated relative resistance to be around 64% when disturbances were larger [23]. Conversely, the chosen value of r involves a conjecture about perceptual sensitivity that is specific to this model and has not been measured directly.

Tietäväinen et al. [18] obtained a value of τ around 300 ms, while other methods of assessment have produced estimates including 125 ms [24] and 200 ms [21]. Direct measurements of ankle response, however, found response to start at 30 ms with maximal displacement around 120 ms [25]. If feedback delay is too long, then intermittent periods of acceleration due to gravity or muscle feedback will be consequently prolonged even as the state enters unstable regions of the phase space. One result is overcompensation for error, in which the fast oscillations found in sway are more amplified than would be the case with shorter delays. Alternatively, if the value of a is too high, then delayed feedback may bypass the unstable manifold and activate at inappropriate locations in the phase space, potentially amplifying slower oscillations over time. Long feedback delays can therefore contribute to instability, sway amplification, and higher risk of falling, but the exact kinds of error are determined by the joint behavior of several parameters, including a, r, and disturbance magnitude σ [1].

Estimation

The above equations represent a Stochastic Delay Differential Equation (SDDE). The Kalman-Bucy filter provides minimum-variance unbiased estimates of the state of a stochastic process when both measurement and process noise are present and can be used to estimate the parameters of continuous-time differential equations from noisy data [26]. However, two challenges arise when estimating the parameters of an SDDE, including the lag interval τ and the lagged position and velocity coefficients, P and D. First, interpolation of the lagged states must be used to allow a continuous domain of possible values for τ. Second, backward extrapolation must be used to estimate the unmeasured interval of lagged states preceding initial state x₀.

Last, we address problems that occur when the discrete switching conditions are toggled between measured instances. For most intervals between measures, the dynamics are linear and the prediction is exact, but state predictions that traverse the condition thresholds will systematically introduce bias to the linear dynamics unless the correct ratio of active and inactive dynamics within each traversal is estimated. We detail an algorithm to resolve this bias by adjusting the prediction according to each of the possible threshold-traversal scenarios.

Optimal filtering

The state-space equation for the time-lagged IPC system is given as

Matrix A contains the parameters of the passive, instantaneous forces, while A_τ contains the conditional parameters of active feedback. When the conditions given in Eq 3 evaluate to false, A_τ = 0.

The measurement matrices simply attribute the observed COM to the state position with estimated origin μ and measurement error variance ϵ:

The complete algebra for the prediction and correction steps of Kalman Filtering is excluded, as its derivation can be found in many resources [27] and remains largely unchanged for this model. However, the key difference in this case is that the prediction step is altered to include the delayed term. Using the following matrix discretizations,

For stationary series with large number of observations, P_t ≈ P_∞. For convenience, we use P_t−1 as an approximation to P_t−τ. Note that Eq 8 does not work if K = 1, making A singular. However, small, numerically viable deviations from K = 1 will not substantially impact solution topology. Point singularitieswill also not impede derivative-free optimization methods.

Estimation of feedback delay

Linear interpolation To obtain estimates of the state at time lags that do not fall on measurement instances, we use linear interpolation of the state:

Backward extrapolation of initial values By introducing an initial value parameter for acceleration, we can estimate a quadratic extrapolation backward from t₀ to t₀ − τ, allowing the influence of lagged states and switching conditions to be respected within the first λ iterations of filtering:

Constrained interpolation of dynamic switching points

To avoid bias due to missing transitional information between measures that straddle the threshold of the conditions given by Eq 3, we explicitly detect each case, interpolate the state falling on the condition threshold, and predict its traversal in two steps. For convenience, take the shortened terms u and v as the delayed states leading up to, and away from the condition threshold:

Conditions for switching off:

If v 2 + v ˙ 2 ≤ r 2, then the lagged state is traversing into the slack radius, and the interpolated point is

In all other cases in which (16) holds true, u is traversing the axis at u = 0.

Conditions for switching on: For cases where the delayed feedback is switching on, the roles of u and v are simply traded. The interpolated point is calculated identically under each set of conditions analogous to those for switching off.

Prediction for threshold traversal: The time for u to reach the switching threshold, Δt⁻, and the time to reach the next observation after the threshold, Δt⁺, can be calculated from the interpolated state at the threshold and its neighboring states, u and v:

In the first step, A, A_τ, and Q are discretized for the interval Δt⁻, and the prediction is given as:

In the second step, A, A_τ, and Q are discretized for the interval Δt⁺, and the prediction is computed from time t + Δt⁻ as:

For either step, A τ d = 0, depending on whether the active feedback is switching on or off.

Optimization

The toggling of active feedback is not a smooth process and results in discontinuities in the space of a cost function for fitting the model, though these are greatly mitigated by the interpolation measures described above. The complexity of the model nonetheless gives rise to multiple local solutions, and attempts to find optimal parameters using local methods such as gradient descent and Nelder-Mead reliably fail. Instead, we recommend using a method of global, derivative-free optimization such as Differential Evolution (DE) [28]. The optimization parameters that we chose are listed below.

• Strategy: DE / rand / 1 / bin with per-vector-dither

• Iterations = 15000

• Population size = 30

• Crossover Probability (CR) = .95

• F = .15

• Weighting of successful members (c) = 0

• Step tolerance: 500

• Relative tolerance: 1e-10

We chose a high crossover probability (CR) due to high dependence between parameters of the model and used simulations to confirm reasonable convergence given the chosen population size, iterations, and F value. DE does not require initial values for parameter estimation, but instead populates a region within explicit bounds. The bounds used here for simulation and data analysis are given in Table 2. Parameter bounds were generally restricted to potentially stable and theoretically meaningful ranges, such as for K, P, r, and τ. Theoretical interpretations of parameters B and D were less certain and were therefore allowed to vary beyond boundaries imposed under any particular physiological definition. τ was constrained to the extremes of the empirical distribution of neural delay given the results from Peterka [21]. Otherwise, bounds were made extreme enough to capture all reasonable possibilities without unnecessarily slowing convergence.

Software

All analyses used R statistical programming environment [29]. Differential Evolution was provided by the R package DEoptim [30]. The IPC model was implemented in C++ using R packages Rcpp [31] and RcppArmadillo [32], and compiled to the open-source R package IPCmodel. The package includes the following functions:

• ipcModel(): C++ Kalman Filter with delayed terms and switching conditions that returns a -2Log-likelihood value for optimization.

• ipcSimulate(): C++ numerical integrator that generates simulated data for the IPC model.

• ipcMultiGroup(): R wrapper for ipcModel() that incorporates physical constants, parameter algebras, and enables the estimation of both within and between-series parameters.

• kalmanIntegrate(): C++ helper function that accepts continuous-time state-space matrices and returns discretized matrices for Kalman-Bucy filtering.

Simulations

Two simulations were used: the first to check model specification, and the second to evaluate the accuracy and precision of parameter estimates. The first simulation used noiseless (i.e. deterministic trajectories) with perfect measurement to check for systematic bias due to the estimation strategy. The second simulation used data simulated to include both process and measurement noise according to the possible properties of real data recorded by a force plate. Solutions for both deterministic and noiseless simulations were generated in linearized steps of size 10⁻⁵s then downsampled according to the design of each simulation. This procedure ensured both numerical accuracy of solutions and simulated real world mapping of analogue processes to discrete measurements. The statistical properties of simulated series were expected to be invariant to downsampling due to the fractal property of continuous random walks (i.e., Wiener processes) where Δt ∼ N(0, Δt).

Parameter sets

Six sets of simulated parameters were defined to test the model’s estimation capability over a variety of possible behaviors and are shown in Table 3. The first set replicates the simulated data for Model 4 by Asai et al. [1] and is named accordingly. The second and third sets (“Low Noise” and “High Noise”) respectively decrease and increase the variance of process noise to examine its effect on other parameters. The fourth and fifth(“Active Control” and “Passive Control”) sets respectively increase and decrease the ratio of active to passive control, representing different plausible configurations for stability. The sixth set (“Rambling and Trembling”) represents a stationary random-walk series that diverges markedly from the underlying theory but is nonetheless a stable and plausible configuration.

Sim 1: Noiseless series

To test for improper model specification and systematic sources of bias, noiseless series were generated to span 20s, with a step size of 10⁻⁵s, then downsampled to an observation every 0.01s and again to every 0.1s. The noiseless series used in the first simulation are shown in Fig 1. Only one series per set and per downsample rate was used, as there were no sources of sampling error. To ensure that estimates converged to a high degree of precision, 3000 iterations of optimization were used.

Table 4 contains the parameter estimates for these simulations, with sampling rate shown to the left. Only the velocity coefficients B and D exhibited substantial bias all throughout, and the “High Active” set incurred the greatest bias over nearly all parameters. Most parameter estimates given 100Hz sampling were exact to at least 3-5 decimal places. Reducing sampling resolution by a factor of ten increased biases to parameters B and D by a factor of ten to fifteen, but much less so for K and P. The nonlinear parameters a, r, and τ exhibited the least bias for all sets.

The small biases to K, B, P, and D most likely occur as a result of the approximate, linear interpolation methods and inability to account for process noise before t₀ in the quadratic backward extrapolation. Biases may be further mitigated using polynomial interpolation of the lagged state. However, the exact accuracy of the estimated τ indicates that bias from linear interpolation is probably trivial in this case.

A second source of bias may be the limits of numerical precision. When no noise is present in the system, the state only occupies a small area of the phase space where certain values of B and D may have nearly unobservable effects on the solution. We show later that relatively unbiased estimates of B and D can indeed be obtained as a function of the other parameters, including the process noise variance σ.

Sim 2 Estimation from noisy data

To test the precision and accuracy of IPC parameter estimates given the dimensions and expected structure of the data from Santos et al. [20], one-hundred individuals were simulated for each parameter set in Table 3, with examples series shown in Fig 2. Each individual consisted of three trials, and each trial consisted of a 60s series downsampled to 100 Hz. The same parameters were estimated for all three trials, making for a total of 18,000 observations per individual model.

Figs 3 through 8 show the sampling variation and bias for each parameter set. Boxplots are grouped by common axis scale. Table 5 gives the means and standard deviations of each parameter for each set.

Variance and baises across all parameters were highly interdependent. Estimates of both process noise (σ) and measurement error were precise and close to their true values, indicating successful filtering of the state. The precision of active and passive and active control parameters depended on the true values of parameters and resulting behavior of the process. For the Asai et al. replication and the sets with low and high process noise, most control parameters had only small bias and high precision, while others were less reliable under particular conditions. The greatest apparent contrast may be the insensitivity radius r, which was not estimable for the Rambling and Trembling set in which its true value was small, and much less reliable in the increased noise set where its value matched Asai et al.

The B and D parameters were the least reliable, and are possibly empirically unidentified without a sufficiently high process noise variance. This is evident from the increased noise set (Fig 6) and the Rambling and Trembling set (Fig 4). The active control set (Fig 8) also showed successful estimation of the B parameter, and improvements in estimating D over the the Asai et al. set, low noise, and passive control.

From the variation in results across sets, it can be inferred that a parameter can only be estimated reliably when the state occurs for a sufficient amount of time in the portions of the phase space for which that parameter has an influence. For instance, the insensitivity radius will not be estimable if the state tends to bypass it entirely. This may be due to a large variance of process noise, or for large values of B that distort the saddle shape of the passive attractor space, causing an orbital path that never intersects the origin. Likewise B and D cannot be estimated reliably if the process does not frequently visit the extremes of the phase space where their influence is most apparent.

Empirical under-identification of some parameters is not necessarily problematic for the others, and does not imply the unreliable parameters should be fixed to some value or excluded. Two solutions to empirical under-identification are to increase the length and resolution of the sample to increase the chances of observing informative behavior, and perhaps to introduce small interventions or disturbances such that subjects express the full range of relevant dynamic behaviors.

Standing apart from the other parameters is the feedback delay τ. Despite perfect accuracy in the noiseless case, it tended to bias downwards when estimated from noisy data. It is unclear from our simulations why the bias occurs and whether it accounts for bias to other parameters. However, the estimates were not generally boundary cases, and the sampling variability was small. If the bias is consistent, the delay parameter should still be comparable between persons, with the caveat that the estimate is understated by 20-40ms.

Data analysis

The IPC model was fit to empirical postural control data to 1) estimate the multivariate distributions of each parameter, 2) test for expected effects from age and visual feedback, 3) test the consistency of parameters within-person, 4) compare the proposed model to simpler alternatives. COM data were obtained from the data set published for public use by Santos et al. [20] and included 49 individuals at 100Hz for 60 seconds per trial. Three trials were conducted with eyes open, and three with eyes closed. Only trials tested with a rigid floor were used for our analyses. Height and weight were provided for each individual and included as the constants h and m in the model, scaled to units of meters and kilograms respectively. Height was scaled by 0.51, the approximate ratio of vertical COM to total height in upright standing (calculated from Table 1, p.7 of [23]). By visual inspection of the sample, it was found that the first and last several seconds of many series contained large, sudden changes in position likely relating to movement during the initiation and termination of the trial period. To ensure that only the stationary dynamics of interest were modeled, 500 occasions were trimmed from the beginning and end of each series, leaving 5000 occasions or 50 seconds of data per trial, and 30,000 measurements in total per individual.

Models

Three models were fit to each of three trials per individual to examine the statistical significance of the parameters involved in intermittent activation and delayed feedback. The models included, in descending order of complexity, the complete intermittent stochastic delay differential equation (ISDDE),

All models included trial-specific initial conditions x_0,i and x ˙ 0 , i and sway origins μ_i for i ∈ [1, 2, 3]. The ISDDE and SDDE both included trial-specific estimation of x ¨ 0 , i for backward extrapolation. All models included measurement error variance σ_ϵ. Parameter boundaries, shown in Table 6 reflected both theoretical and analytic roles of each parameter. For example, B could not be less than zero in the ISDDE because it is conjectured to represent ankle stiffness, and stability is required to come from values of P and D in the given domains. In the SDE, stable solutions must rely on only instantaneous feedback with coefficients K and B. In the absence of other theoretical mechanisms, the same physiological interpretations of K and B could not be assumed and thus the same theoretical constraints were not applied.

Multiple regression was used to test the association between each parameter, visual feedback, and age, accounting for height and mass as covariates. Pearson correlation was used to estimate the correlation between parameter estimates during trials with eyes open and trials with eyes closed. Maximum likelihood estimation was used to fit each model, assuming the multivariate normality of measurement and process noise.

The estimated means μ ^, standard deviations σ ^, and medians of each estimated parameter across all trials × participants × visual feedback conditions, are given in Table 7. The estimated individual-level intraclass correlations (ρ_ICC), effect sizes, and p-values for age and visual feedback are also given for each model. Measurement error estimates were generally small (σ_ϵ < 1e − 3) and were omitted from the tables. Minor, trial-specific “nuisance” parameters including sway origins and initial values were also omitted. Mean sway origin was estimated to be 0.217, with a standard deviation of 0.11 and a median of 0.226.

The estimated parameters of the ISDDE fell within the expected domains. Several parameters of the ISDDE had outliers that substantially inflated estimates of their standard deviations. Trimmed standard deviations in which the highest 15 values were excluded are given in parentheses in Table 7. The marginal distributions of each parameter with these trimmed means and standard deviations are shown in Fig 9. B, P, D, and r in particular were skewed upward by outliers but otherwise had relatively precise distributions about their medians, with similar precision to those of the SDDE. K had consistent values around .91 to .93 in all three models. B was close to zero for most series but skewed upward by outliers as high as 80. In the SDE, B was allowed to take negative values but had a mean around 17. All values of B in the SDE were positive and greater than zero, with a minimum of .82. a was generally high, representing active control over 75-85% of the phase space. Similarly, r was 50% smaller on average than values used in previous studies. τ had a median of 284 ms and was distributed between 200 to 400 ms. If the bias found in simulations is consistent and proportional, then the true median delay was closer to 240 ms. The SDDE estimated much longer delays on average at 470-490 ms but much lower values of D. Process noise standard deviation σ_w estimates were distributed identically between models.

No significant effects of visual feedback were observed in the parameters of any of the three models. The lowest p-values were for a(p = .069) and σ_w(p = .086). Both the SDDE and SDE showed effects on σ with p < .05, the alpha level before adjusting for the 17 tests in total.

In the SDDE, both passive ankle stiffness K and active control coefficient P were shown to significantly increase with age. The effects were detected given the adjusted alpha level, with p < .0029. Both B and σ_w in the SDE showed significant trends with age as well. Other non-significant effects with p < .05 were K and σ_w in both the ISDDE and SDDE, B in the ISDDE, and τ in the SDDE. Effect estimates of σ_w were consistent across models.

Overall, parameters tended to be more consistent within person for the simpler models. The highest intra-class correlation for all parameters in all models was K, with extreme reliability (ρ = .999) in the SDDE. σ generally correlated around .5 for each model. The ISDDE had the least consistent parameters with intraclass correlations near zero for P, D, a, and r. The SDDE and SDE intraclass correlations were moderate to high for all except feedback delay, τ, which was near zero.

Akaike’s Information Criterion (AIC) [33], as - 2 ln ( L ^ ) + 2 k where k is the number of estimated parameters, was used to compare overall model fit for every individual. For each trial, the model with the lowest value of the AIC was selected as the best fitting option. In total, the ISDDE was selected for 227 trials, SDDE for 62, and SDE for 0. No significant associations were found between model selections over trials within person or by visual feedback condition.

Discussion

Simulation

Simulation studies were used to determine whether the parameters of the intermittent activation feedback control model proposed by Asai et al. [1] can be estimated using a Kalman Filtering-based framework with delayed proportional and derivative terms and discrete activation thresholds. The results of the simulations show that the parameters of the model can be estimated with relatively low bias and high precision if the behaviors for which they are influential are sufficiently expressed in the data (i.e., empirically identified). Every parameter of the model was successfully recovered in at least one of the parameter configurations tested, though no single configuration of parameters resulted in a completely unbiased set. The set of results shown in Fig 8 comes close, with downward bias only to the active derivative controller. We can also see by comparing Fig 6 to Figs 5 and 3 that an increased variance of process noise allowed the identification of the B and D parameters, but with large standard errors. The derivative coefficients were likely biased and unreliable when the trajectories did not frequent the extremes of position and velocity where the directional effects of derivative terms could be distinguished from other sources. Fig 4 shows that with process noise of a standard deviation much greater than the insensitivity radius and a weak attraction to point equilibria (K ≈ 1), the state is prone to drifting away from the origin where it will rarely traverse the insensitivity radius or switching boundary. If the data can be optimally explained without the use of the switching parameters, then they are said to be empirically unidentified. For this reason, both a and r do not contribute crucial information and converge to precise solutions when the data are optimally described by other parameter values characteristic of rambling and trembling. However, estimates of the derivative controllers B and D in that case were unbiased.

Across configurations, some iterations of model fitting resulted in negative values of B. In continual PID controllers, this would result in amplification and instability over time. In the ISDDE and SDDE, the stability of the system given a value of B depends on the corresponding values of P, D, and a, as the instantaneous proportional and derivative terms do not control the complete periodic behavior on their own. Negative values of B will promote further instability in the already unstable manifold of instantaneous feedback but will be counteracted when the state reaches the stable manifold determined by active feedback. It is informative that solutions occasionally involved negative values of B that breach its theoretical interpretation as joint friction. Solutions that did better identify B and D only did so when their effects were much larger than physically plausible, a priori values of joint friction. In simpler PID cases, estimates of damping tend to be far less reliable than, for instance, the proportional coefficients, so for these reasons together it may be inadvisable to rely on postural sway data and estimation approaches to specifically determine joint friction. Similar concerns may be directed toward the active feedback damping D, though the prior ISDDE literature does not assert as specific of a definition nor necessary theoretical boundaries.

Estimates of both process noise and measurement error were very close to their true values in every case, with only small upward bias proportional to the magnitude of the estimate for certain parameter sets (Figs 3 and 4). The standard deviation of measurement error that we chose to simulate was σ = .01 cm, twenty times the error of the force plate used by Santos et al. [20] to obtain the data. The success of estimation despite greatly exaggerated sensor noise demonstrates the reliability of Kalman filtering and adequate technical specification of the model, and relieves researchers from the need to choose a preliminary noise reduction step such as spectral filtering. Instead, using the raw data and including measurement error in the model avoids removing fine-grained details of the signal represented in the domain of high frequencies typically suppressed by low-pass filtering.

The feedback delay, τ, was unbiased in noiseless simulations, and consistently biased downward in noisy simulations. It is not clear what causes the bias, but it did not appear to consistently induce bias in other parameters that depended on the correct lag interval, such as the active proportional and derivative controllers.

The results of our simulation demonstrate that the proposed method of direct, statistical estimation by Kalman filter can recover the complete set of parameters for the model. Previous estimation methods only attempted to estimate five of the eight structural parameters [17, 18]. Among those attempted, the D parameter did not converge to its true value in simulation nor to a reliable, unimodal distribution in the empirical study. Despite this setback, no discussion was given of the role of empirical identification in determining D or other parameters, whereas we have demonstrated that the precision of estimates depends on their true values and interdependence. Additionally, the accuracy of their results rests on assumed values of K, B, and r. Due to the high degree of parameter dependence in univariate models such as this, error in one parameter is expected to propagate to other parameters in a compensatory manner. It is therefore preferable to jointly estimate all uncertain model parameters when possible.

The prior studies also did not account for measurement error. We determined that additional sources of sensor noise could be filtered simultaneously with estimation of the dynamic structure. If additive noise is Gaussian, then no preprocessing steps such as spectral filtering or downsampling should be needed and the risk of obscuring important, fine-grained topological features is greatly mitigated.

Computationally, the use of global optimization to maximize the likelihood function provided an efficient alternative to Bayesian MCMC methods as no data simulation procedures, prior distributions, or posterior sampling were required. An additional, unexplored benefit of maximum likelihood in this case is estimation of standard errors directly from the likelihood function. Because the model includes discrete thresholds, the likelihood function was stochastic and non-differentiable. This prevented the use of the Hessian matrix to calculate precision. However, future work may explore methods of smoothly approximating the marginal likelihood function, for instance by fitting splines to likelihood values retained from the optimization procedure.

Experimental data

The results of analyzing the empirical COM data show that the nonlinear mechanisms of feedback activation led to significant improvement in model fit over the simpler SDDE and SDE (i.e., delayed and instantaneous PID) models. It cannot be determined from statistical model comparisons alone whether the results validate the model-generating theory of posture control. To that end, we must compare the parameter estimates to their theoretical priors.

Overall, the distributions of parameters showed a feasible correspondence to the domains expected given the theory. K was consistently close to the 91% relative resistance found by Loram and Lakie [22] for all of the models tested, here showing resistance to 92% of the total gravitational toppling torque on average. Conversely, in the ISDDE and SDDE, B most often converged to zero and was not likely to play a critical role in the model behavior. Perhaps coincidentally, the mean of B was near its proposed value of 4 Nms/rad. It is possible that statistical power at the individual level was insufficient to identify small effects due to B, and the expected value would be recovered if it were estimated across the total data set. Active feedback was generally weaker than hypothesized but still sufficient for stability. Estimates of P were closer to .1 than the proposed .25 [1], likely due in part to the greater resistance to toppling forces from values of K closer to the high end of their theoretical distribution. D played a large role in the dynamics of active control and resulted in non-negligible damping in many individuals. Values of a and r reflected greater control sensitivity than expected. a values around 75% to 80% assign a larger share of the phase space to active feedback, while smaller values of r indicate less tolerance to falling at the origin of sway. The mean estimate of a was higher than found by Tietäväinen et al. [18], which reported a control space closer to 64% in accordance with the analysis by Asai et al. [1]. We found a nearly identical distribution of the feedback delay, τ, to Tietäväinen et al. [18], ranging from 200 to 400 ms with a mean around 300 ms. Estimates of σ_w were an order of magnitude smaller than expected by Asai et al. [1], and about half of those found by Tietäväinen et al. [18].

A graphical vignette of these results is provided in Fig 10, which shows six raw data series with their respective intermittent activation conditions estimated by the model. The horizontal axis is the tilt angle and the vertical axis is the tilt angle’s velocity. The shaded region represents behavior where P and D are equal to zero. In the unshaded region, all parameters are active with their non-zero values. Fig 10a shows two cases that resemble the theorized structure with combinations of stable and unstable manifolds in nearly equal proportion. In Fig 10b, the values of B and a are sufficiently large to minimize the influence of the unstable manifold. The result is behavior that closely resembles harmonic oscillation around a single equilibrium. The opposite trend is shown in Fig 10c, where the unstable manifold is not influential, but a high ratio of the derivative coefficients to proportional coefficients results in continual suppression of velocity. This pattern results in wandering oscillations without clear equilibria.

The optimal solution for the SDE model had a much larger, positive value for B than the other models while maintaining a theoretically plausible value of K less than 1. Furthermore, all values of B in the SDE were positive, as we might expect given that negative values would result in instability. When a linear system is strongly overdamped (in this case, a high value of B in the SDE, or either B or D in the SDDE) with relatively weak proportional feedback, as the SDE, then it exhibits non-equilibrium Langevin dynamics. These dynamics have conventionally described the random walk of a large molecule due to its collisions with a many smaller molecules in a solvent. The resulting trajectories can appear locally stationary by chance and exhibit short intervals of oscillation. Previous studies have modeled posture control in the context of Langevin dynamics [34–37]. Our simplest model of COM movement, the SDE, resembled a model of COP proposed by Bosek et al. [34] that describes trajectories as a second order SDE with no proportional feedback and a large derivative coefficient B. Fig 11 shows how the theoretical model and Langevin dynamics differ markedly in their mechanistic parameterization and observed phase portraits, yet they share many notable features. In both, high-frequency oscillations move gradually across the sample space in a “rambling” pattern. By chance, the Langevin equation in Fig 11b can result in concentrated oscillations around a few apparent equilibria, but no equilibrium mechanism is present in the model. The parsimony of generating these patterns with only three parameters poses a challenge to the specificity of evidence for the theoretical ISDDE. Visual inspection of the complete results showed that trials ranged between the two extremes of theoretical misspecification, from harmonic oscillation to Langevin dynamics. The expected topology involves a mix of features from both, sometimes showing adherence to the principles of feedback switching with occasional deviations into Langevin-type random walk.

Regardless of the true form of the underlying process, we might expect that if the parameters represent underlying physiological mechanisms, they should exhibit some degree of trait-like stability within-person. The intraclass correlations in Table 7 show that the nonlinear switching parameters were generally unreliable within-person. The correlations for the remaining parameters increase as the model is simplified to the SDDE and SDE. The higher consistency of the simpler models’ parameters does not necessarily imply that they are more “real” than those of the ISDDE. It is expected that reliance on fewer parameters to explain the variance of sway results in fewer competing configurations of those parameters. Any consistency of topological features within-person will be reflected in similarity of the model solutions. The lack of consistency in the more complex ISDDE is, however, a challenge to the trait-like stability and actuality of its parameters.

Though no specific connections between visual feedback and the theoretical mechanisms of control were hypothesized for the present study, we expected one or more parameters to be significantly influenced over trials in which eyes were closed in correspondence with previously observed effects on summary statistics. By modeling center of pressure variation with Langevin dynamics, Bosek et al. [34] found that the process noise distribution was influenced by visual feedback. The same finding was replicated with further connections to Parkinson’s disease [35]. Vieira et al. [38] found associations of visual feedback with stabilogram measures of sway. All three models models tested here had lower p-values for σ_w than for other parameters, suggesting that effects may be discernible given a larger sample or improvement in model specification.

Age has been previously associated with more general metrics of sway, such as path length [39], frequency band [38, 40] and mean velocity [40], though findings vary and few effects have been consistently reproduced in COP and COM data. Significant effects of age were observed in the present results, including ankle stiffness and active feedback force in the SDDE and process noise and ankle viscosity in the SDE. Interpretation of effects on the SDE is more difficult because the parameters of the SDE do not correspond to specific explanatory mechanisms in this case. The consistent positive associations of all models with noise magnitude σ_w with age may be linked to previously observed associations of stabilogram-based diffusion metrics with age [41]. The significant associations of ankle stiffness and active proportional feedback in the SDDE found here may reflect previously observed increases in stiffness and damping with age estimated from a simpler PID model [42].

Finally, the model concerns an abstract notion of body tilt angle, though there are many ways to represent this using the full kinematic data. For simplicity and consistency with past studies, we chose to represent tilt angle by the COM. Preliminary tests using alternative measures included COP and the average angle of both ankle joints. The results were found to differ markedly from both our current results and those previously obtained with the COM, but a complete comparison of alternative measures is too complex to discuss here. We leave detailed examination of this question with regard to the feasibility of this model to future study.

Conclusions

We designed and implemented an Extended Kalman Filter-based estimation model of intermittent, delayed feedback control in postural sway and demonstrated that for a variety of stable configurations, parameters can be recovered accurately given adequate empirical identification. Application of the model to experimental data resulted in distributions of the parameters the correspond well to previous findings and suggest that physiologically informative and clinically useful attributes of human balance may be extracted directly from COM data. While the model replicates previous findings, the conjectured parameters of feedback activation were not reliable within-person or strongly associated with visual feedback and age. Further comparisons with alternative mechanistic theories and model parameterizations are warranted. Beyond postural control, the model stands as a framework for estimating parameters of stochastic delay differential equation models controlled by discrete activation thresholds.