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

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.


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 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: where I is the rotational inertia, m is the body mass (kg), g is gravity (� 9.81m/s 2 ), and h is the height of the COM. T includes all the terms representing mechanisms of resistance to the angular toppling force. w t is a Gaussian, independent and identically distributed random variable accounting for stochastic variation in acceleration, with standard deviation σ. The total passive forces may be written as mgh(1 − K)θ t , as K is the percentage of the gravitational acceleration counteracted by passive resistance. While a certain definition of B is not given, its effects are non-trivial and an interpretation may be taken from the common use of the secondorder damped oscillator equation, in which the velocity coefficient represents negative feedback due to friction. In this case, it may be regarded as a measure of ankle and knee joint friction.
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 a ¼ sinðapÞ cosðapÞ , we change the upper and lower estimation bounds from [−1, 1] 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 0 .
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 where Q is the process noise covariance matrix, H is the measurement matrix, μ is the estimated origin about which the COM oscillates, and R is the covariance matrix of measurement error. The contemporaneous and lagged state vectors are ; and the state transition matrices are 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, we can then provide the prediction equations for the state mean and covariance as follows: For stationary series with large number of observations, P t � P 1 . 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: λ is the conversion of the time delay to the number of measured occasions comprising that interval. The ceiling and floor functions thus give valid measurement indices and are used to give a combination of measurements falling to either side of λ, weighted proportionally. If τ = 0, then the second term of Eq 13 can be neglected. Backward extrapolation of initial values By introducing an initial value parameter for acceleration, we can estimate a quadratic extrapolation backward from t 0 to t 0 − τ, 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: Where s is the seconds constant, such that v; u; _ v, and _ u are measured in radians. For use later, we note here that the slope between the two points is

Conditions for switching off:
then A τ is switching off. If this holds true, then the following conditions further apply: and !½ðv > 0 and _ v < 0 and u < 0 and _ u < 0Þ then the lagged state is traversing the line _ x ¼ ax outside of the slack radius and not traversing u = 0. The interpolated point ðû;_ uÞ falls on the line, and is calculated aŝ 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.
then A τ is switching on. If this holds true, then the following conditions further apply: then the lagged state is traversing the line _ x ¼ ax outside of the slack radius and not traversing u = 0, and the interpolated point is calculated as Eq (18). If u 2 þ _ u 2 � r 2 , then the lagged state is traversing the slack radius from within, and the interpolated point is calculated with Eq (19). In all other cases in which Eq (21) holds true, u is traversing the axis at u = 0 and the interpolated point is calculated as Eq (20).
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:x In the second step, A, A τ , and Q are discretized for the interval Δt + , and the prediction is computed from time t + Δt − as:x For either step, A d t ¼ 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. 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.

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 −5 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 −5 s, then downsampled to an observation every Direct estimation of parameters of an intermittent feedback model of postural sway 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 0 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

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 a stochastic delay differential equation (SDDE) with delayed feedback but no intermittent switching conditions, All models included trial-specific initial conditions x 0,i and _ x 0;i and sway origins μ i for i 2 [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, Direct estimation of parameters of an intermittent feedback model of postural sway 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.   Table 3, and parameter descriptions are given in Table 1. The estimated meansm, 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.

Par
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Þ þ 2k 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.

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][35][36][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 traitlike 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 traitlike 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]. Phase portraits and time series from two models generated from the same vector of noise (scaled by σ w ). In (a) data were simulated from the ISDDE model with the theoretical priors. In (b) data were simulated from a 3-parameter SDE. The large ratio of derivative (B) to proportional (K) force results in non-equilibrium Langevin dynamics that may exhibit similar features to the ISDDE for limited periods of time. Characteristic features distinguish (a) from (b), such as sharp changes in velocity localized to quadrants II and IV, slow change in velocity associated with quadrants I and III, and higher density at two spatial equilibria. https://doi.org/10.1371/journal.pone.0222664.g011 Direct estimation of parameters of an intermittent feedback model of postural sway 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.