Abstract

We developed a theory of human stance control that predicted (1) how subjects re-weight their utilization of proprioceptive and graviceptive orientation information in experiments where eyes closed stance was perturbed by surface-tilt stimuli with different amplitudes, (2) the experimentally observed increase in body sway variability (i.e. the “remnant” body sway that could not be attributed to the stimulus) with increasing surface-tilt amplitude, (3) neural controller feedback gains that determine the amount of corrective torque generated in relation to sensory cues signaling body orientation, and (4) the magnitude and structure of spontaneous body sway. Responses to surface-tilt perturbations with different amplitudes were interpreted using a feedback control model to determine control parameters and changes in these parameters with stimulus amplitude. Different combinations of internal sensory and/or motor noise sources were added to the model to identify the properties of noise sources that were able to account for the experimental remnant sway characteristics. Various behavioral criteria were investigated to determine if optimization of these criteria could predict the identified model parameters and amplitude-dependent parameter changes. Robust findings were that remnant sway characteristics were best predicted by models that included both sensory and motor noise, the graviceptive noise magnitude was about ten times larger than the proprioceptive noise, and noise sources with signal-dependent properties provided better explanations of remnant sway. Overall results indicate that humans dynamically weight sensory system contributions to stance control and tune their corrective responses to minimize the energetic effects of sensory noise and external stimuli.

Introduction

The control system for human bipedal upright stance involves the generation of an appropriately calibrated corrective torque based on body-sway motion detected primarily by vestibular, visual, and proprioceptive sensory systems (Horak and Macpherson 1996). Because body motions are small in this task, one expects that the signal-to-noise ratios of sensory signals are poor. Based on previous studies of sensory integration (Ernst and Banks 2002), one might hypothesize that the nervous system generates corrective torque based on an optimal estimate of body orientation derived from a weighted combination of the noisy sensory cues (Fig. 1, left side Sensory Integration component). From Bayesian estimation theory, optimal sensory weights can be determined that provide a sensory representation that is a maximum likelihood estimate of a physical variable, S, if the variances of the different sensory signals are known. An appropriately weighted combination of the noisy sensory signals and will provide an optimal estimate, , of the physical variable with lower variance than either of the individual sensory signals. For the Sensory Integration component shown in Fig. 1 and with the assumptions that the sensory noise sources are independent, Gaussian, and the Bayesian prior is uniform, then the optimal sensory weights can be calculated from the variances of the sensory signals.

Schematic depiction of the stance control system and factors potentially influencing its control. Considering the sensory integration component of the system in isolation (system in the dashed box), Bayesian estimation theory shows that optimal sensory...

However, the application of this estimation method to the stance control system is more complex. The variability of any particular sensory signal depends not only on the inherent noise properties of the individual sensory systems, but also is confounded by the system operating with feedback (Fig. 1, right side Motor/Biomechanics combined with the left side Sensory Integration). This confound occurs because variability that might initially arise from one particular sensory system contributes to the production of a time varying body sway, which in turn is sensed by the other sensory systems and, therefore, contributes to the variability of the other sensory signals. Additionally, the overall variability of sensory signals is influenced by the possible contribution of a motor noise component (Harris and Wolpert 1998) and by the filtering properties provided by neural feedback gains, muscle (activation) dynamics, and body dynamics (van der Kooij et al. 1999). That is, with feedback and the addition of motor noise, the assumption that the variability in different sensory signals is due to independent noise sources is not satisfied. Therefore, the relatively simple calculation of a maximal likelihood estimate of shown in Fig. 1 for the Sensory Integration component of the system operating without feedback would not directly apply to a system operating with feedback although the general concept of sensory weighting may still be valid.

The presence of a continuously applied external perturbation adds further complications by contributing a stimulus-evoked component to the overall variability of sensory signals in addition to inherent variability in each sensory system’s signal. For example, stance on a continuously rotating surface directly contributes to increased variability of proprioceptive signals which encode body motion relative to the surface. This proprioceptive orientation signal generates body sway that is sensed by the graviceptive system (mainly vestibular) and the visual system (if vision is available). Therefore, variability of neural signals from graviceptive and visual systems is indirectly increased by the presence of an external surface tilt stimulus. Furthermore, the variability of an overall internal sensory orientation estimate that consists of a weighted summation of orientation signals from individual sensory systems will depend on the nervous system’s choice of sensory weights for each sensory channel.

Despite these complications, it is reasonable to hypothesize that the nervous system is able to account for the intrinsic variability in sensory and motor systems, the feedback nature of stance control, and the variability contributed by external perturbations in order to generate body sway behavior that is in some sense optimal. The goal of the current study is to determine whether an optimal estimation theory can be developed that accounts for these various complications and that explains a wide variety of features of experimental results from a previous study that characterized stimulus-evoked sway in humans (Peterka 2002) and additional features related to body-sway variability derived from a re-analysis of the previous data.

Specifically, we sought to develop an optimal theory of stance control that provides a parsimonious explanation for three key experimental findings. These are (1) the systematic changes in stance control dynamics that are dominated by a decrease in body sway response sensitivity to stimuli of increasing amplitude, (2) the systematic increase in “remnant” body sway variability observed with increasing stimulus amplitude where remnant sway variability is defined as the variability that is not directly attributable to the applied stimulus, and (3) the fact that the corrective torque generated in relation to body motion appears to be specifically selected to achieve some goal that involves more than just the need to maintain stability.

Materials and methods

We used a multi-stage approach to investigate the hypothesis that a model described below (Fig. 2) and associated optimizations based on this model can account for the experimentally observed sensory re-weighting, the properties of the neural controller, and the patterns of remnant sway. The analysis stages include (Stage 1) estimation of stance control system parameters and sensory weights, (Stage 2) exploration of alternative internal noise sources (sensory only, motor only, and combinations of sensory and motor sources) and various functional forms of these noise sources to determine which noise source best accounts for the experimentally observed remnant sway, (Stage 3) given a particular noise source, identification of a behavioral criterion that predicts the observed sensory re-weighting, (Stage 4) given a particular noise source, identification of a behavioral criterion that predicts the observed neural controller parameters, and (Stage 5) demonstration that the identified noise source accurately predicts spontaneous sway behavior.

Model of human stance control in which sensory information from proprioceptive and graviceptive systems is weighted (wp and wg) to provide an estimate of body-in-space sway angle, bsest. However, bsest is potentially biased from true body-in-space orientation...

Experimental study

Experimental data used in the current study were part of an extensive experimental data set previously used to characterized the human stance control system under a variety of conditions by applying perturbations that evoked sagittal plane body sway by tilting a visual surround viewed by a subject and/or tilting the support surface (ss) upon which a subject stood (Peterka 2002).

Data from 8 healthy subjects (age range 24–46 year) who participated in this previous study were used in the current study. All subjects gave written informed consent to a protocol approved by the Institutional Review Board of Oregon Health & Science University.

The current study made use of only a portion of the previous data for the condition where subjects stood with eyes closed with tilt perturbations applied to the ss at five different amplitudes (0.5, 1, 2, 4, and 8° peak-to-peak) on different trials. Test trials in the prior study were presented in randomized order. A wide-bandwidth pseudorandom tilt perturbation based on a mathematical integration of a step-like waveform derived from a pseudorandom ternary sequence (PRTS, (Davies 1970; Peterka 2002)) was used to evoke the sagittal-plane center-of-mass (CoM) body-in-space sway angle (bs). The PRTS was created using a 5-stage shift register operation with feedback to create a repeating sequence of numbers with each number having a value of 0, 1, or 2. The feedback was chosen to produce a maximal-length sequence prior to repeating (242 sequence length). This sequence was transformed into a sampled stimulus time series by mapping each number to a value representing the ss stimulus velocity (0➔0, 1➔+v, and 2➔ -v°/s) and holding each velocity value for Δt = 0.25 s. This ss velocity waveform was mathematically integrated and scaled in amplitude to give a position waveform that was sampled at 100 Hz and used to command the servo-motor controlling the angular tilt of the ss. The PRTS had a period of 60.5 s. Eight complete cycles were presented for the 0.5° stimulus and six cycles were presented for the other amplitudes. A stimulus based on a maximal-length PRTS has the property that only odd-numbered spectral components have nonzero energy and the non-zero components of the ss velocity waveform have approximately constant amplitude up to a frequency of ~1/(3Δt) = 1.3 Hz.

For the data used in the current study, all tests were performed with subjects using a backboard assembly that constrained their body mechanics to be that of a single-link inverted pendulum. Results from the previous study demonstrated that use of the backboard did not alter postural dynamics in this particular experiment (Fig. 7 in (Peterka 2002)). The mass, moment of inertia, and CoM height of the backboard assembly were taken into consideration in all of our modeling efforts.

Results of the Stage 4 analysis. The shaded areas indicate the derivative (combination of neural controller kd and intrinsic damping bi) and proportional (combination of neural controller kp and intrinsic stiffness ki) gains for which the stance control...

Experimental data analysis

Body sway and stimulus waveforms were decomposed using a discrete Fourier transform (DFT) into frequency-component parts at frequencies (f) ranging from 0.017–1.3 Hz. The frequency response functions (FRFs) from ss to bs at the five stimulus amplitudes (indicated by subscript a), Hss2bs,a(f), were calculated at each excited frequency by dividing the mean (across stimulus cycles) DFT values derived from bs by the mean DFT values derived from ss. The complex numbers representing Hss2bs,a were expressed as gain and phase values as a function of the frequency. FRFs were calculated for each subject and then averaged across all subjects since the FRFs were previously shown to have similar form despite variations in subject’s anthropometric measures (Peterka 2002). See Appendix A for details of equations used for the calculation of experimental FRFs.

The single-sided power spectral density (PSD) of remnant body sway at each stimulus amplitude, Pbsr,a(f), was calculated using methods that exploit the periodic nature of the PRTS (van der Kooij and de Vlugt 2007), see Appendix A. Remnant PSDs were calculated for each subject and then averaged across subjects to reduce the variance of the final estimates of the Pbsr,a(f) spectra. At the odd-harmonic frequencies that correspond to frequencies excited by the PRTS stimulus, the remnant PSD values represent the variance of the body sway that was not attributable to the stimulus-evoked sway. At the even-harmonic frequencies that were not excited by the stimulus, the remnant PSD values represent the variability in body sway at these frequencies.

Stance control model

We developed a feedback control model to aid in the interpretation of the experimental results (Fig. 2). The model is a slightly modified version compared to the model originally proposed in (Peterka 2002). The differences include: 1) the PID (Proportional, Integral, Derivative) controller was replaced by a PD controller, and torque feedback was added that better accounts for the low frequency gain declines and phase advances seen in experimental FRFs (Peterka 2003; Cenciarini and Peterka 2006); 2) hypothetical sensory and motor noise sources were added; and 3) muscle activation dynamics were included to represent the filtering action of the conversion from motor control signals to muscle force. For the eyes closed condition, the model includes feedback from proprioception that signals body sway relative to the surface, from graviception (mainly vestibular origin) that signals body sway relative to earth vertical, and from somatosensory receptors that sense the torque acting on the body to control stance. The signals from proprioception and graviception are assumed to provide accurate, wide-bandwidth (i.e. no dynamics) orientation estimates. As such, these signals do not represent the dynamic characteristics of primary sensors, but rather represent processed sensory representations of body motion that are likely available for stance control (Angelaki et al. 1999; Merfeld et al. 1999; Casabona et al. 2004; Bosco et al. 2005). Sensory integration was represented as a weighted combination of the proprioceptive and graviceptive contributions with weight factors wp and wg, respectively. These weight factors quantify the relative contribution from these two sensory systems such that wp + wg = 1. These combined sensory signals form an internal body-in-space orientation estimate, bsest, which is compared to an internal reference of bs = 0 (not shown in Fig. 2) to produce an error signal that contributes to the generation of corrective torque. We previously found that low-frequency FRF characteristics could be explained well if we assumed there was also a feedback contribution from a low-pass filtered, torque-related sensory signal (Peterka 2003; Cenciarini and Peterka 2006). This torque feedback can be thought of as providing an additional, but in this case a slowly time-varying reference signal, bsref, that produces an error signal and subsequent corrective torque that drives the steady-state value of bs to an orientation that requires very little corrective torque. For example during a sustained surface tilt, an error signal based only on bsest would result in a sustained body tilt away from vertical with greater body tilt for larger values of wp. However when the torque feedback loop is included, the body will be driven slowly back toward an upright vertical orientation. Thus, the system can maintain an upright stance orientation on a tilted surface while still taking advantage of a weighted combination of proprioceptive and graviceptive cues.

The signal representing the combined proprioceptive, graviceptive, and torque sensory information drives a “neural controller” that generates corrective torque, tc, that includes a position-related component (kp gain factor) and a velocity-related component (kd gain factor). A time delay is included in the control loop that represents the combined delay due to various subsystems (sensory transduction, sensory processing, neural transmission, muscle activation). The corrective torque is applied to the ankle joint of an inverted pendulum body to produce body sway. Appendix B includes the equations of inverted-pendulum body dynamics, muscle-tendon, activation, and neural control dynamics and also the transfer functions that were used in the different stages of analysis.

In this paper, “transfer functions” refer to the Laplace transform domain representations of the model-derived differential equations that characterize the dynamic relationship between the various input variables in the model (e.g., ss signal, noise signals, etc.) to output signals (most often bs). The frequency domain representation of a transfer function is easily obtained, and can be displayed as gain and phase functions versus frequency. The model-derived transfer functions can be compared, when useful, with experimentally derived transfer functions which we refer to as “frequency response functions” to distinguish between model and experimental results.

Stage 1 analysis: model parameters from fits to FRFs

The goal of the first stage was to investigate whether the FRFs from ss to bs could be described by the stance control model, and whether the systematic change in response dynamics could be captured by only a change in the sensory weights. The stance control model shown in Fig 2 was fitted to the five FRFs to determine a parameter vector containing intrinsic joint stiffness (ki) and damping (bi), the neural controller proportional (kp) and derivative (kd) gains, a neural time delay (τd), the gain (kt) and time constant (τt) of the low-pass filter in the torque feedback loop, and proprioceptive and graviceptive weights (wp and wg) for each of the five ss amplitudes, a. The parameter vector was found by the minimization of a loss function (J1) that was summed across the FRFs and model transfer functions from the five stimulus amplitudes (indicated by subscript a), and the sum across frequency of the squared normalized difference between the experimental FRFs, Hss2bs,a (f), and the model transfer functions of surface stimulus to body sway, Hss2bs,a(f,θ1):

The differences between the FRFs and the model transfer functions were 1) normalized by the FRF magnitude to insure that differences in frequency regions where FRFs have lower gain values were just as important as differences between high gain values and 2) weighted by the inverse of the frequency to prevent over-fitting of the higher frequency data where, on a logarithmic scale, there were more data values.

Stage 2 analysis: noise models to explain remnant sway

The goal of the second stage of analysis was to identify intrinsic noise models that could account for the body sway remnant PSDs, Pbsr,a(f) estimated from experimental data. Fits included the levels and shapes of sensory and motor noise spectra. For the noise spectra we initially considered both white noise filtered by a first order low-pass filter and noise with a spectral distribution of 1/fα noise, i.e. white (α = 0), pink (α = 1), or Brownian noise (α = 2). The results of the low-pass filtered white noise models in all cases provided poorer fits to remnant PSD than the 1/fα models, and were therefore disregarded in the results presented in this study.

The level of sensory and motor noise was a combination of baseline motor noise level and, dependent on the noise model, a signal dependent component expressed by a noise-to-signal ratio (NSR) scalar r. The NSR was either defined with respect to the mean-square value of the signal or to a frequency-dependent PSD.

The noise model parameter vector was determined by minimization of a loss function (J2) that is the sum of the squared normalized differences between the estimated and model body sway remnant across all stimulus amplitudes:

At each stimulus amplitude, the PSD of the model-predicted body sway remnant is the sum of and , which represent the body sway remnant PSDs produced by sensory and motor noise sources, respectively. The difference between the experimentally measured remnant PSDs and model-predicted remnant PSDs was normalized by the experimental remnant PSD and weighted by the inverse of the frequency.

Sensory noise models

In the case of sensory noise, the PSD of the model-predicted body sway remnant is shaped by the transfer function from the sensory noise sources to body sway, and by the spectral properties of sensory noise sources:

where Pvt, Pvp, and Pvg are the PSDs of torque, proprioceptive, and graviceptive noise, respectively. The subscripts of the transfer functions denote the input and output of the transfer functions that are given in Appendix B. For example Hvt2bs is the transfer function from torque signal noise to body sway. We implemented different models of the sensor noise spectra (Pvt, Pvp, Pvg).

Signal-independent sensory noise

In the first sensory noise model (S1) the noise spectra of the different sensory systems have the same spectral shape and only differ in level. The spectra for S1 are given by:

where the level of the noise is denoted by the scalar n and the t, p, and g subscripts of n denote the torque, proprioceptive, and graviceptive sensory signals, respectively. In sensory noise model S1, Ps(f) is a PSD that is the same for all sensory noise signals, and has a spectral shape characterized by 1/fα noise:

For the S1 model, the noise-model parameter vector θ2 included the parameters nt, np, ng, and αs to be identified by minimization of the J2 loss function.

Signal-dependent sensor noise

We also considered cases where the noise spectra scaled with the average signal level (mean-square value) or with the spectrum of the signal evoked by the external surface tilt stimulus.

The spectra for sensor noise model S2 is the sum of the S1 spectrum and an additional part that scales with the sensor spectrum induced by the external stimulus. Specifically, the torque, proprioceptive, and graviceptive sensor noise spectra are given by:

where the rs’s are scalar noise-to-signal ratios that multiply the PSDs Pt(f), Pp(f), and Pg(f) that denote the stimulus-dependent portion of the overall sensory noise spectra (see Appendix B). For the S2 noise model, the noise-to-signal ratio rs was the same for every sensory system. An alternative model was also investigated that allowed different rs values for each sensory system, but optimization results suggested that a model of this form was over-parameterized.

For sensory noise model S3, the shape of the spectrum resembles 1/fα noise as in S1 but the noise level also scales with the average mean-square value of the sensory signal induced by the external stimulus. The three sensor noise spectra are given by:

where E{()2} denotes the mean-square value of the sensory signals evoked by the external stimulus. These mean-square values can be calculated from the corresponding spectrum (see Appendix B). For S3 noise model, the noise-to-signal ratio rs was the same for every sensory system. As with the S2 model, a model of the form of S3 that allowed different rs values for each sensory system appeared to be over-parameterized.

Motor noise models

In the case of motor noise, the spectral shape of the body-sway remnant PSD attributable to motor noise (Pbsr_m) is determined by the transfer function from the motor noise source to body sway (Hvm2bs) and the spectral shape of the motor noise (Pvm):

For the motor noise spectrum we considered different models. In all models the motor noise spectrum was signal dependent. The different motor noise models included a baseline level (nm) and a part that scaled linearly (rm) with either the mean-square value or the spectrum (Ptc) of the corrective torque signal induced by the external stimulus.

For the first motor noise model (M1), the spectrum resembled 1/fα noise that scaled with the mean-square value of the corrective torque induced by the external stimulus:

In the second motor noise model (M2), the spectrum was the sum of a signal-independent component with a 1/fα spectral shape and a signal-dependent part that scaled with the spectrum of the corrective torque induced by the external stimulus. For M2 the motor noise spectrum is given by:

The third motor noise model is a combination of M1 and M2.

In M3 the noise level scaled both with the mean-square value of the signal and with the spectrum at the excited frequencies. This was done to account for the characteristics of the experimental body-sway remnant PSD that showed an increase in noise at the non-excited frequencies with increasing stimulus amplitude but an even larger increase at the excited frequencies.

The goal of the Stage 3 and 4 analysis stages was to identify a behavioral optimization criterion that could predict the identified control-model parameters identified in the Stage 1 analysis. More specifically, we aimed to predict the control-model parameters that could be regulated by the nervous system to achieve some desired goal: these parameters are the sensory weights (Stage 3 analysis) and the neural controller feedback gains (Stage 4 analysis). To this end, a loss function related to behavioral criteria was minimized:

We consider six different behavioral criteria that included four “performance” criteria and two “effort” criteria. We defined performance criteria as kinematic criteria related to body-sway motion. In the case of minimization of a performance criterion, is given by:

is the model-predicted mean-square value (determined by calculating the mathematical integral of the various PSDs) of either body-sway angle (n = 0), velocity (n = 1), acceleration (n = 2), or jerk (n = 3). As in any linear system, the output spectrum (in this case body sway) is the sum of all input spectra multiplied by the squared absolute transfer function of the corresponding input signals (sensory noise, motor noise, and the external tilt stimulus) to the output (Bendat and Piersol 2000). The term 2πfn is a n-order differential operator in the frequency domain applied to the spectrum of body sway, in order to obtain spectra of body sway velocity, acceleration, or jerk.

We defined an effort criterion as a kinetic criterion related to the generation of corrective torque. In the case of minimization of an effort criterion, is given by:

For n = 0, is the mean-square value of the corrective torque and for n = 1, is the mean-square value of the derivative of corrective torque (referred to as torque change). Similar to the calculation, the mean square value of torque (n = 0) or torque change (n = 1) was obtained by integrating the model-predicted PSD of corrective torque or torque change. The PSD of corrective torque or torque change was the sum of all input spectra multiplied with the squared absolute transfer function of the corresponding input signals to the output.

For the Stage 3 analysis the goal was to find which of the six behavioral criteria could predict the change in sensory weights due to increasing ss stimulus amplitude. For each of the six behavioral criteria, the parameters, , identified by the Stage 3 analysis were the sets of sensory weights (wp’s and by definition the wg’s since wg = 1–wp) that minimized the model-predicted mean-square value ( or ). In the calculation of and , the control-system parameters, , identified in the Stage 1 analysis (except for the sensory weights) and the noise-model parameters, , identified in the Stage 2 analysis were held fixed.

For the Stage 4 analysis the goal was to find which of the six behavioral criteria could predict the specific values of neural controller feedback gains. For each of the six behavioral criteria, the parameters, , identified by the Stage 4 analysis were the neural controller positional (kp) and derivative (kd) gains that minimized the model-predicted mean-square value ( or ). In the calculation of and , the control-system parameters, , identified in the Stage 1 analysis (except for the neural controller gains) and the noise-model parameters, , identified in the Stage 2 analysis were held fixed.

Stage 5 analysis: prediction of spontaneous sway characteristics

The goal of the Stage 5 analysis was to investigate whether the identified stance control model along with the identified internal noise sources was able to predict the structure of spontaneous sway. Spontaneous sway is the body sway during quiet stance when there is an absence of external perturbations. A random walk analysis of spontaneous body sway reveals a structure that can be visualized by a stabilogram diffusion plot (Collins and De Luca 1993). In such a plot the mean squared difference between two samples are given as a function of the difference in time between two corresponding samples. Two regions have been identified in stabilogram diffusion plots. In the short term region, corresponding to time differences less than ~1 s, body sway can be characterized as a persistent (unstable) random walk. In the long-term region body sway can be characterized as an anti-persistent (stable) walk.

The model identified in Stage 1 was simulated in Matlab Simulink with input noise disturbances defined by the simplest sensory noise model (S1) from the Stage 2 analysis. The S1 noise model was used because this model resulted in reasonably good fits to the experimental body sway remnant PSDs (see Section 3) and was easy to simulate. Although the results from the Stage 2 analysis showed that some other sensory noise models and combinations of sensory and motor noise models resulted in even better fits to the remnant PSDs, we reasoned that if S1 gave reasonable predictions of spontaneous behavior, the more complex noise models would give similar results since the differences between S1 and the more complex noise models were small. For the fitted sensory weights we used the condition closest to quiet stance condition, i.e. the condition with the smallest support surface amplitude. The sensory noise signal was generated from a white noise signal that was transformed to the frequency domain using a DFT. This spectrum was shaped according to the identified S1 sensory noise, and transformed back to the time domain. The model-derived stabilogram diffusion function was calculated from the time series of body center-of-pressure (CoP) displacements that were calculated from the simulated body sway. A mean experimental stabilogram diffusion function was also calculated by averaging across the eight diffusion functions of individual subjects calculated from 360 s recordings of anterior-posterior CoP data obtained in eyes closed conditions. Curve fits to the model-derived and experimental stabilogram diffusion functions provided estimates of the short and long term diffusion coefficients, scaling coefficients, and critical coordinates that were compared to previously reported values (Collins and De Luca 1995).

Results

Stage 1 analysis

The experimentally determined gain and phase values of the mean FRFs showed systematic changes as a function of the amplitude of the surface tilt stimulus (Fig 3(a), points connected by dotted lines). The FRF gain values at each individual frequency generally decreased with increasing stimulus amplitude such that the sets of gain values corresponding to the different stimulus amplitudes maintained very similar shapes across the frequency range. The FRF phase data from different stimulus amplitudes had similar values at frequencies of about 0.1 Hz and below but showed some divergence at higher stimulus frequencies with the largest phase lag for the lowest stimulus amplitude and the least phase lag for the highest stimulus amplitude.

The model transfer function curve fits to the experimental FRFs in the Stage 1 analysis provided a good description of the experimental data that captured the major features of the gain and phase changes with frequency and as a function of stimulus amplitude (Fig. 3(a), solid curves). Note that in the Stage 1 analysis, the internal sensory and motor noise characteristics were not taken into account since they do not affect the estimated FRF or the model transfer function between the external surface tilt stimulus and the body sway response. Most of the identified parameters values (given in Fig. 3 legend) were similar to those previously reported for the particular set of parameters (kp, kd, ki, bi, τd, and wp) that were included in both the current model (Fig. 2) and the previous model (Peterka 2002). However, the time delay parameter, τd = 97 ms, from the current analysis was noticeably shorter than the time delay previously reported (about 170 ms). This is explained by the fact that the previous model did not include muscle-activation dynamics. Specifically, both muscle-activation dynamics and a time delay impart a phase lag at higher frequencies. Inclusion of muscle-activation dynamics accounts for some of the phase lag observed at higher frequencies and the time delay that is needed to account for the remainder of the phase lag is reduced.

The Stage 1 analysis demonstrated that a good description of the experimental FRFs could be obtained even though the analysis imposed a strong constraint that, with the exception of wp and wg = 1–wp , no model parameters could vary as a function of stimulus amplitude. The success of the model description in accounting for the experimental FRFs indicates that the Fig. 2 model provides a parsimonious interpretation that attributes all of the amplitude-related changes in FRFs to a sensory re-weighting phenomenon whereby the proprioceptive (wp) and corresponding graviceptive (wg = 1–wp) contributions change with stimulus amplitude. In particular, the identified wp decreased with increasing stimulus amplitude (Fig. 3(b)). Again, this pattern of wp change with stimulus amplitude is similar to that previously identified (Peterka 2002) using a slightly different model and allowing all model parameters in individual fits to vary across subjects and stimulus amplitudes.

Stage 2 analysis

The goal of the Stage 2 analysis was to determine whether internal sensory and/or motor noise sources could be identified that were able to account for the experimentally observed remnant power spectra of body CoM sway variability and were compatible with the previously identified model parameters identified in the Stage 1 analysis.

As shown in Fig. 4(a), both the remnant sway and the stimulus-evoked sway increased with increasing stimulus amplitude. The frequency distributions of the remnant sway are represented by the family of power spectra shown in Fig. 4(b). The remnant PSD values at each individual frequency generally increased with increasing stimulus amplitude and the sets of PSD values corresponding to the different stimulus amplitudes maintained similar shapes across the frequency range when plotted on log-log scales. The PSDs generally decreased in magnitude with increasing frequency. For frequencies below about 0.3 Hz, the PSDs decreased approximately in proportion to f-1, and for frequencies above about 0.7 Hz there was a steeper decrease that was approximately proportional to f-4.

(a) The body-sway variability can be decomposed into a part that is evoked by the stimulus and into the remnant sway that is not directly attributable to the applied stimulus (Appendix A). Both the experimental stimulus-evoked and remnant body...

A detailed feature of the remnant PSDs is the saw-tooth shape of the PSDs. This saw-tooth structure is most evident in the PSDs at lower frequencies and higher stimulus amplitudes where the PSD values at even-frequency harmonics tended to be lower in magnitude than the adjacent odd-frequency harmonics. The PRTS stimulus only had power at odd frequencies and these were the frequencies that tended to have the greater remnant power. However, it is important to understand that there is no a priori expectation that the remnant PSD will have greater amplitudes at frequencies where there is stimulus energy as compared to frequencies where there is no stimulus energy. Therefore, the presence of the saw-tooth pattern is indicative of some coupling between the stimulus and the remnant body sway. A likely source of this coupling is a contribution from internal noise sources that have a signal-dependent component, meaning that the magnitude of the noise depends, to some extent, on the magnitude of sensory and/or motor signals within the stance control system.

The Stage 2 analysis considered potential noise sources that included only sensory noise (proprioceptive, graviceptive, and force/torque noise), only motor noise, or combinations of sensory and motor noise with various plausible functional forms (see Methods). Both signal-independent and signal-dependent noise sources were investigated. The sensory noise sources were added to the sensory system signals in the Fig. 2 model prior to sensory integration by linear summation of the proprioceptive and graviceptive signals, and prior to the low-pass filtering associated with the processing of the force/torque signal. Thus, these sensory noise sources represent variability due to peripheral encoding as well as variability due to central processing of peripheral signals within a given sensory system. The motor noise source contributed to the overall corrective torque and represents the variability associated with generating the desired corrective torque by muscle activation.

The Stage 2 analysis used various combinations of sensory and/or motor noise models in combination with the stance control model parameters identified in the Stage 1 analysis to predict the remnant PSDs shown in Fig. 4(b). For the various types and combinations of noise models investigated, the parameters of the noise models were adjusted to minimize a loss function to obtain optimal fits to the family of remnant sway PSDs.

Table 1 summarizes the minimum value of the loss function for different sensory and motor noise models, and combinations of these models. The uniformly large loss function values for all of the motor-only noise models (M1-M3) indicate that none of these models provided adequate fits to the remnant PSDs. All of the sensory-only noise models (S1-S3) provided much lower loss function values than the motor-only noise models, but several noise models that included both sensory and motor noise (S2M1, S1M3, S2M3, and S3M3) yielded still lower and similar loss function values of about 100. Both motor and sensory noise models that included signal-dependent noise terms were able to capture the saw-tooth structure of the remnant PSDs.

The results of four representative fits of noise models to the remnant PSDs are shown in Fig. 5. One includes only motor noise (M3 top left), one includes only sensory noise (S1 bottom left), and two include combinations of motor and sensory noise (S2M1, S1M3 right column). All of the noise models predicted remnant PSDs that increased in magnitude with increasing stimulus amplitude. However, there were large differences between noise models in the extent to which the various models could account for the stimulus-amplitude dependency, the variation across the full frequency range, and the detailed saw-tooth structure of the experimental remnant PSDs. The remnant PSDs predicted from the motor-only noise models were similar to the experimental PSDs only in the mid-frequency range of about 0.1 to 0.3 Hz, and even in this range the prediction for the highest stimulus amplitude was uniformly lower than the experimental data. At both lower and higher frequencies, the motor-only models predicted much lower PSD values across all stimulus amplitudes than experimentally measured.

Results of the Stage 2 analysis. Examples of different noise model predictions of remnant power spectra for models with sensory only, motor only, or combinations of sensory and motor noise. The model predictions were derived from noise model fits to the...

The simplest sensory-only noise model (S1) provided a much better fit (Fig. 5) to the remnant PSDs than any of the motor-only noise models. The PSD fits from this model were quite close to the experimental data for all stimulus amplitudes over frequencies ranging from about 0.06 to 1.3 Hz (the highest frequency evaluated). Only the model predictions for the lowest frequencies and the highest stimulus amplitudes appeared to be inconsistent with the experimental data. Additionally, because this simple model did not include a signal-dependent noise component, the model prediction was unable to account for the saw-tooth structure of the PSDs. The sensory-only models that included signal-dependent noise (S2, S3; plots not shown) were both able to account for the saw-tooth PSD structure, but there were only modest reductions in their loss functions compared to the S1 model.

The four models with combined sensory and motor noise and with the lowest loss functions (S2M1, S1M3, S2M3, and S3M3) all provided similar fits to the remnant PSDs (two of these fits are shown in Fig. 5) and all included signal-dependent terms that account for the saw-tooth nature of the PSDs. Although the quality of these fits and their associated loss functions are similar, the parameters of these four fits (Table 1) demonstrate that it is difficult to reach a strong conclusion as to whether the signal-dependent noise should be attributed to a motor or sensory source. In S2M1 and S2M3 the sensor noise is highly signal-dependent in comparison to the motor noise. But in S1M3 and S3M3 the motor noise is highly signal-dependent in comparison to the sensory noise.

All of the sensory noise models included parameters that defined the general magnitude of noise associated with the sensory systems. The noise magnitude estimates for proprioceptive noise, np, and graviceptive noise, ng, were similar across all noise models. A consistent finding was that the noise magnitude in the graviceptive system was about 11 times greater than the noise in the proprioceptive system (mean ng/np = 11.1 ± 0.83 SD for the four models with the lowest loss function values). This disparity in noise between these two sensory systems accounts for the large increase in remnant sway as the stance control system shifts toward greater use of graviceptive cues (increasing wg) and decreasing use of proprioceptive cues (decreasing wp) with increasing stimulus amplitude.

Estimates of torque noise magnitude, nt, were consistent across the three models that included only sensory noise (S1-S3, Table 1). However, the torque noise estimates became quite inconsistent across noise models when both sensory and motor noise components were included. This inconsistency is attributable to the fact that motor noise and the low-pass filtered torque noise both accounted similarly for the enhanced amplitude of remnant sway at low frequencies such that removal of both motor and torque noise from the model resulted in remnant sway predictions that severely underestimated the actual remnant sway at low frequencies (results not shown). The interaction between motor and torque noise was evident in that the noise model fits that assigned large values to nt also assigned small values to nm, the parameter representing the fixed portion of motor noise. Conversely, fits with small nt values had large nm values. Thus, the Stage 2 analysis indicated that it was necessary to include some fixed internal noise source that contributed to low frequency remnant sway, but the analysis was not able to determine whether this noise source had a sensory or a motor origin.

Stage 3 analysis

The goal of the stage 3 analysis was to determine whether the pattern of wp and wg changes with stimulus amplitude identified in the Stage 1 analysis could be predicted. The hypothesis was that the stance control system selected wp and wg values to minimize a behavioral criterion. Six different behavioral criteria were investigated that included the minimization of the mean-squared value of body CoM sway angular position, sway angular velocity, sway angular acceleration, sway angular jerk, corrective torque, and the rate-of-change of corrective torque. Minimizations of these behavioral criteria were used to predict the wp values at the 5 stimulus amplitudes. The stance control model parameters (except for the wp and wg parameters) were held fixed at values previously identified in the Stage 1 analysis, and the internal noise properties were set to those that provided good fits to the remnant sway PSDs in the Stage 2 analysis.

The wp Stage 3 predictions overlayed with the wp values from the Stage 1 analysis are shown in Fig. 6 for the six behavioral criteria. The wp predictions were essentially identical for all of the sensory/motor noise models from the Stage 2 analysis that gave similar predictions of the remnant sway PSDs (cost function values of about 100). The particular wp prediction shown in Fig. 6 is from the S1M3 noise model. The wp predictions for the sensory only noise models were also nearly identical to the results shown in Fig. 6 (wg not shown since wg= 1 - wp).

All of the behavioral criteria predicted decreasing values of wp with increasing stimulus amplitude. This is intuitively expected when the proprioceptive noise is much smaller than the graviceptive noise. At low stimulus amplitudes there is little stimulus-evoked sway, so the overall magnitude of any of the behavioral measures, which is due to both stimulus-evoked and internal noise components, is smallest when the system is relying primarily on the lower-noise proprioceptive system. However, with increasing stimulus amplitude, the stimulus-evoked sway would become quite large if wp remained the same. The overall sway and all of the behavioral criteria can be reduced if the stance control system shifts toward increased reliance on the graviceptive cues (even though the graviceptive system is much noisier than proprioception), thus reducing the responsiveness of the system to surface perturbations.

While all of the behavioral criteria produced the correct trend of decreasing wp with increasing stimulus amplitude, the sway velocity behavioral criterion provided wp predictions that were closest to the wp values identified in the Stage 1 analysis. Sway acceleration, sway jerk, and rate-of-change in torque criteria all overestimated wp. The sway position and torque criteria both provide good wp predictions at the lowest stimulus amplitudes (0.5° and 1°) but underestimated wp at the larger stimulus amplitudes.

In contrast to wp predictions based on sensory or combined sensory/motor noise models, models that included only motor noise were completely unable to predict wp changes with stimulus amplitude (Fig. 6). Specifically, models with only motor noise predicted that wp should be zero independent of the stimulus amplitude. This result is intuitively expected because, with wp = 0, the surface-tilt stimulus evokes the smallest possible body sway for any given stimulus amplitude.

Stage 4 analysis

The goal of the stage 4 analysis was to determine whether the neural controller kp and kd parameters identified in the Stage 1 analysis could be predicted. The hypothesis was that the stance control system selected specific kp and kd values to minimize a behavioral criterion. The same six behavioral criteria investigated in the Stage 3 analysis were also applied to the Stage 4 analysis. Minimizations of these behavioral criteria were used to predict the kp and kd values at the five stimulus amplitudes. The stance control model parameters (except for the kp and kd parameters) were held fixed at values previously identified in the Stage 1 analysis, and the internal noise properties were set to those that provided the best fits to the remnant sway PSDs in the Stage 2 analysis.

For all six behavioral criteria, the Stage 4 analysis produced kp and kd values that did not noticeably change with stimulus amplitude. This is consistent with the Stage 1 results which showed that it was not necessary to change the values of kp and kd parameters as a function of stimulus amplitude in order to account for the amplitude-dependent changes in FRFs (Fig. 3).

Minimizations of the six different behavioral criteria led to six different sets of kp and kd parameters. Each of these six sets of kp and kd values, combined with the passive stiffness (ki) and damping (bi) factors from the Stage 1 analysis, are plotted in Fig. 7 along with the kp + ki and kd + bi values obtained from the Stage 1 analysis. Additionally, the plot shows the three different nested regions that correspond to the ranges of kp + ki and kd + bi values that are compatible with stability of the stance control model system for three different time delay values. The largest region is for the shortest time delay of 100 ms which is close to the time delay of 97 ms identified in the Stage 1 analysis. The stability regions associated with increasing time delay values become progressively smaller (regions for 150 ms and 200 ms time delay values are shown). For time delays of above about 340 ms, no kp and kd values can be selected that provide stability.

Minimizing the mean squared value of torque provided the closest prediction of the kp and kd values identified in the Stage 1 analysis (Fig. 7). The sway velocity minimization provided kp and kd predictions that were almost as close to the Stage 1 results as the torque prediction. The predictions from the other four behavioral criteria were noticeably more distant from the Stage 1 results.

Stage 5 analysis

The goal of the stage 5 analysis was to demonstrate that the stance control and noise models, which were characterized entirely by analysis of stimulus evoked sway, could also predict spontaneous sway behavior. Quiet standing center-of-pressure (CoP) displacement was simulated using the model parameters identified in the Stage 1 and 2 analyses. The results of the simulation were used to calculate a stabilogram diffusion function of the CoP movements in the anterior-posterior direction, and to measure the short and long term diffusion coefficients and the critical point coordinates that describe its characteristic shape (Fig. 8(a)). The model-derived stabilogram diffusion plot showed the typical shape described previously (Collins and De Luca 1993). That is, the function increased linearly with increasing time interval until the critical time point at Δtc, and then continued to increase linearly, but with a lower slope, at long time intervals. The results of the simulation were compared with the average across-subject stabilogram diffusion function (Fig. 8(b)) calculated from 363 s recordings of anterior-posterior CoP data obtained in eyes closed conditions. The diffusion coefficients and critical point coordinates of the model simulations and experiments resembled one another. The values of these parameters were in the range of those estimated from 30 s records of anterior-posterior CoP in eyes closed conditions (Collins and De Luca 1995), except for Ds and ΔXc, which were larger in our study for both the simulated and the experimental results. These differences may be due to the sensitivity of RMS measures of spontaneous sway to sample durations whereby the RMS sway increases with sample duration (Carpenter et al. 2001). Since our experimental data records were more than ten times longer than the duration of earlier studies, the sway magnitude, as reflected by the Ds and ΔXc parameters, was also larger.

Results of the Stage 5 analysis. Shown are the stabilogram diffusion functions, SDFs, and their diffusion coefficients (short term, Ds, and long term, Dl) derived from linear fits (dotted lines) to the SDFs, and critical coordinates (critical time, Δ...

Discussion

Sensory re-weighting mechanisms have been shown to contribute to human stance control (Kiemel et al. 2002; Peterka 2002) as well as to other sensorimotor tasks (Safstrom and Edin 2004; Mugge et al. 2009). Previous investigations indicate that an optimal weighting of information from multiple noisy sensory sources can be used to formulate an overall estimate that has lower variance than an estimate formed from any of the individual sensory sources, and that re-weighting can occur to compensate for conditions that degrade the accuracy of a particular sensory source (Ernst and Banks 2002). However in a sensorimotor feedback control system such as the stance control system, the inherent noise of sensory systems is only one factor influencing an overall body orientation estimate and therefore the variability of body sway. Other factors include the contributions of motor noise and the variability caused by the continuous application of a noise-like external perturbation. Optimal estimation and control theory predicts that an improved orientation estimate can be obtained by down-weighting of the particular system that is perturbed (van der Kooij et al. 2001). However in the stance control system, down-weighting of one sensory system must be compensated by up-weighting of another to insure there is sufficient torque generated to resist gravity (Cenciarini and Peterka 2006). Depending on the sources of internal noise, the re-weighting that compensates for an external perturbation should also affect the variability of the remnant body sway, and thus give insight into the internal sources of variability.

We showed that a stance control model (Fig. 2) accounted very well for the experimental response dynamics and remnant sway. This model included two parallel feedback pathways (proprioceptive and graviceptive feedback), an inner torque feedback loop, and muscle-tendon dynamics. A systematic decrease in the proprioceptive weight and corresponding increase in the graviceptive weight was sufficient to account for the decrease in the experimentally measured FRF gain that characterizes the sensitivity of body sway evoked by the surface rotation stimulus.

Both sensory and motor noise sources were investigated as possible sources for body sway variability. We compared different models of sensory and motor noise, and combinations of these models. Noise models that included only motor noise were clearly inferior. Models that included only sensory noise were far superior to the motor-only noise models, but were inferior to all models that included both sensory and motor noise. There was no single model that included both sensory and motor noise that was clearly superior to all others since several models resulted in similar remnant fits. However, there were similarities among the best combined sensory and motor noise models: 1) sensory noise provided the major contribution to remnant body sway, 2) the noise included a signal-dependent component although it was unclear if the signal dependency was best attributed to sensory or motor noise; 3) the baseline level of graviceptive noise was about ten times larger than that of proprioceptive noise; and 4) the sensory noise spectrum resembled 1/fα noise, with α close to one (pink noise). In our approach we attributed the time varying properties of body sway to external stimuli, sensory noise, and motor noise. Possible variations in neural strategy and muscle dynamics were captured by motor noise in our method. Possible variations in sensory dynamics were captured by sensory noise in our method. We note that the estimated sensory and motor noise were not white but pink. Compared to white noise, pink noise has enhanced energy at low frequencies. This enhanced low frequency variability could be thought of as capturing behavior corresponding to slow variations in time that could be due to slow fluctuations in sensory dynamics, neural strategy, or alterations in neuron-pool dynamics.

A previous study of smooth pursuit eye movements provided evidence that sensory noise is a major determiner of behavioral properties related to the initiation of eye movements (Osborne et al. 2005), but during continuous eye tracking the system is apparently dominated by signal-dependent motor noise (Medina and Lisberger 2007). Most literature has focused on systems where motor noise, and particularly signal-dependent motor noise appears to be prominent. The variation in motor unit recruitment thresholds and twitch forces causes muscular force to exhibit monotonically increasing signal-dependent motor noise (Jones et al. 2002). In systems with signal-dependent motor noise, a minimum-variance theory accurately predicts the trajectories of both saccades and goal-directed arm movements, the speed- accuracy trade-off described by Fitt's law (Harris and Wolpert 1998), obstacle avoidance (Hamilton and Wolpert 2002), step tracking wrist movements (Haruno and Wolpert 2005), and the orientation of arm impedance relative to a curved object (Selen et al. 2009).

Compared with other studies, we evaluated more complex models of neural noise. We found that the noise included both a signal-independent (baseline) and signal-dependent components. In most other studies noise consisted of either a baseline level (van der Kooij et al. 1999) or a signal-dependent component (Harris and Wolpert 1998), but not combinations of both. In computational models that considered sensor and/or motor noise, noise was assumed to be lowpass filtered white noise (Harris and Wolpert 1998) or Brownian noise (Liu and Todorov 2007). However, a robust finding in our study was that 1/f sensory noise provided the best fits. 1/f or pink noise is a common feature of scale free networks (Albert and Barabási 2002) and is observed in many biological signals such as heart rate (Saul et al. 1988) and EEG (Freeman 2008). The 1/f noise characteristic is also found in the long range correlations observed in human sway (Duarte and Zatsiorsky 2001).

We found that noise models that explained the structure of remnant sway in stimulus-evoked sway conditions were also able to predict the structure of spontaneous sway as characterized by the characteristic 2-part shape of the predicted stabilogram diffusion function (Fig. 8). Another robust finding was that the baseline level of graviceptive noise was about ten times larger than the baseline level of proprioceptive noise. The large difference in noise levels in these two sensory systems largely explains the increase in remnant sway with increasing stimulus amplitude. Specifically, larger surface motions increase variability in the proprioceptive signal, thus evoking a down-weighting of proprioception and up-weighting of graviception. Although the shift toward greater utilization of graviception decreased the stimulus-evoked sway, the remnant sway was increased due to the larger amplification of graviceptive noise due to the larger graviceptive weighting factor.

One might hypothesize that the nervous system could “know” that a particular type and amplitude of an external perturbation was occurring and therefore would use this knowledge to adjust sensory weights to reduce the influence of the perturbation and to prevent instability. But considering the feedback structure of the stance control system, the presence of internal noise in sensory and motor systems, and variability induced by a continuous external perturbation, one would expect there to be uncertainty in knowledge about the particular perturbation giving rise to the complex, time-varying pattern of sensory signals. However, a mechanism that is able adjust sensory weights to minimize the effect of internal (sensory and motor noise) and external perturbations does not need to know the nature of the external perturbation. Furthermore, because larger amplitude perturbations evoked larger sensory signal variability, the increased variability would drive larger changes in sensory weights (assuming that weight changes are able to reduce overall variability for a particular perturbation). Therefore without assuming that the system has specific knowledge about the perturbation and its amplitude, the weight changes would scale appropriately with perturbation amplitude in order to prevent instability.

All the behavioral criteria investigated predicted a systematic decrease of the proprioceptive weight with increasing surface stimulus amplitude when sensor noise or a combination of both sensory and motor noise were included in the model. These model predictions were generally consistent with the experimental results (Fig. 6). The dominant effect of the sensory noise contribution was demonstrated by observing that when only motor noise was included in the model the optimal solution was to completely ignore proprioception for all amplitudes of support surface motion. Thus the model prediction with only motor noise was not at all consistent with the experimental results and supports the finding that remnant sway results were not predicted by models with only motor noise.

In the case of sensory noise or combinations of sensory and motor noise, the behavioral criterion that best predicted the pattern of sensory re-weighting was minimization of body sway velocity, and thus was not one of the commonly used criteria such as minimization of deviation (Harris and Wolpert 1998; van der Kooij et al. 1999; Todorov 2004), jerk (Hogan 1984; Flash and Hogan 1985), rate-of-torque change (Uno et al. 1989), or control effort (van der Kooij et al. 1999; Carver et al. 2005). No single behavioral criterion we investigated was able to accurately predict the experimentally determined neural controller gains. The experimental neural controller gains were midway between the gains predicted by minimizing body sway velocity and corrective torque. Given that the product of corrective torque and body sway velocity equals mechanical work, it seem likely that mechanical work, and thus metabolic energy, is a likely candidate criterion for minimization during standing. Metabolic energy minimization has also been shown to predict preferred gait parameters (Donelan et al. 2001).

In theory, humans could make predictions about the near future of body orientation and then continuously integrate the delayed sensory measures with predictions to form an estimate of their current body orientation. In visuomotor control it has been shown that the motor output that the brain programmed to start a reaching movement or to correct it midflight was a continuous combination of two streams of information: a stream that predicted the near future of the state of the environment and a stream that provided a delayed measurement of that state (Izawa and Shadmehr 2008). In the past we and others (van der Kooij et al. 2001; Kuo 2005) developed models of human stance control that also included the prediction of sensory and body states to improve the estimation of body orientation. However, in the current model we ignored this prediction of future states since it was not needed to account for our experimental data. Currently there is no experimental evidence that prediction of future states is involved in human stance control, but also no strong counter evidence.

In conclusion, our work supports the theory that during stance humans dynamically weight orientation information from multiple sensory sources and set neural feedback gains to minimize the effect of external perturbations and signal-dependent sensory and motor noise to minimize an energy related criterion. A strong prediction of our results is that sensory re-weighting is dependent on the spectrum of external stimuli since the weighting of sensory clues is a tradeoff between managing the disturbing effects due to internal sensory noise and external stimuli.

Acknowledgements

This work was supported by National Institutes of Health Grant AG-17960 and by the BrainGain Smart Mix Programme of the Netherlands Ministry of Economic Affairs and the Netherlands Ministry of Education, Culture and Science.

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Appendix A

For the estimation of the experimental FRFs relating the ss stimulus to the bs response and the PSDs of remnant body sway we used a non-parametric method that exploits the periodic nature of the PRTS perturbation signal (Pintelon and Schoukens 2001). The input was a periodic PRTS signal with cycle time T = 60.5 s. The input and output signals had a sample time of Ts = 0.01 s and, therefore, each cycle contained Np = T/Ts = 6050 samples. In the experiment we collected six or eight successive cycles (each with length Np). To avoid transient responses, the first cycle was not included in the analysis. Therefore the data set consisted of M = 5 or 7 cycles. By calculating the DFTs of the M data blocks and averaging them we can separate the periodic, stimulus-related component of the bs response from the non-periodic component of the response (i.e., the remnant response). The periodic component is used to calculate the FRF and the non-periodic component is used to calculate the remnant PSD.

The analysis first calculates the DFTs of each cycle’s data:

for l = 1 to M cycles with n = 1 to Np representing the time sample index and k representing the frequency index such that the index k corresponds to the frequency 2π k/T Hz.

The mean DFT across the M cycles is then calculated to determine the periodic component of the stimulus and response:

The experimental FRF between ss and bs is obtained by calculating the ratio:

Note that this ratio was only calculated at odd numbered values of k where there was energy in the PRTS ss stimulus. In the paper we use the frequency variable f as the argument in the experimental FRFs, but it should be understood that f is calculated only at discrete values corresponding to frequencies where there was energy in the PRTS. Furthermore, when comparing experimental FRFs with model transfer functions for the purpose of calculating fit errors and minimizing cost functions, the model transfer functions were evaluated at frequencies corresponding to the frequencies of the FRFs.

The remnant body sway PSDs are calculated by determining the variance of the sampled DFTs:

This calculation is made at both odd and even harmonics of the PRTS and represents the portion of the overall PSD of bs motion that is not attributable to periodic components of body motion evoked by the periodic PRTS stimulus.

In the case of a linear system, Pbsr at non-excited frequencies should only have energy related to the non-periodic component of bs. In the case of a system with nonlinearity, the non-excited frequencies would also include responses to the periodic stimulus that depended on the specific nonlinearity. It is important to understand that the response components due to the nonlinearity would also be periodic. Therefore the remnant PSD would not include periodic components due to the system nonlinearity, but would only include components due to non-periodic noise.

Appendix B

The transfer function dynamics of the different blocks of the human stance control model are defined here using a Laplace transform representation (s = the Laplace variable) of the differential equations of the various portions of the stance control system shown in Fig. 2.

The rigid body dynamics are those of a single-link inverted pendulum:

where J = 81.1 Kg m2 is the moment of inertia around the ankle joint, m = 83.3 Kg the body mass, g = 9.801 m/s2 the gravitational accelerations and h = 0.896 m the distance of the center-of-mass from the ankle joint.

The body is stabilized by a corrective torque that includes contributions from both intrinsic joint and muscle-tendon dynamics and active neurally-mediated control. The muscle/tendon dynamics and passive joint properties produce a visco-elastic torque modeled by a spring (ki) in series with a damper (bi):

The neural controller includes proportional (kp) and derivative (kd) feedback gains and a lumped time delay (τd) that represents all delays in the system:

Prior to being summed with the signal representing torque due to intrinsic muscle/tendon dynamics, the signal from the neural controller is fed through the muscle-activation block that has low-pass filter characteristics:

The input of the neural controller is the difference between an internal set-point and the sum of graviceptive orientation information weighted by wg, proprioceptive information weighted by wp, and information related to the corrective torque, tc, applied at the ankle joint. The sum of wp and wg is always one. The graviceptive signal, g, is a neural representation of the bs angle assumed to be created in the CNS through the fusion of semicircular canal and otolith information (Angelaki et al. 1999). The proprioceptive signal, p, is a neural representation of the ankle joint angle (bs-ss) assumed to be created in the CNS by the fusion of all sensors (probably primarily muscle spindles) capable of contributing to a joint angle signal. The torque sensor signal, t, is a neural representation of the corrective torque applied to the ankle joint. The contribution from force/torque feedback provides an effective way to resist long-lasting disturbances that would otherwise produce a sustained forward or backward lean of the body. This mechanism is implemented by positive feedback of a low-pass filtered version of the torque sensor signal to provide a slowly time-varying set point signal, bsref. The low-pass filter is:

where kt is the feedback gain and τt is the time constant of a low-pass filter (Peterka 2003).

The neural signals representing proprioception, graviception, and torque are assumed to be imperfect in that these signals are noisy. This sensory noise is modeled by adding a noise signal to each of these sensory signals. Because humans are not able to maintain a perfectly constant force (Jones et al. 2002), motor noise is added to the corrective torque in the model.

Predictions of body sway responses to surface tilt stimuli and of sway variability are facilitated by defining the following set of transfer functions that were used in Stages 1-4 of our model-based analysis. First we define the sensitivity function S that will be used to simplify the other transfer functions we derived:

The transfer function that defines how a hypothetical external torque perturbation, tex, applied to the body evokes bs can be derived from the model and is given by:

Note that in the subscript of the preceding and in the following transfer function equations, the ‘2’ denotes ‘to’ and is used to express the direction of the signal flow from the stimulus signal to the response signal. The transfer function of the surface tilt stimulus, ss, to body sway, bs, is given by:

The transfer functions of ss rotation to tc and to the torque sensory signal, t, are given by:

The transfer function of ss rotation to the graviceptive sensory signal, g, is given by:

The transfer function of ss rotation to the proprioceptive sensory signal, p, is given by:

The transfer function of the motor noise, vm, to bs is given by:

The transfer function of proprioceptive signal noise, vp, to bs is given by:

The transfer function of gravioceptive signal noise, vg, to bs is given by:

The transfer function of torque signal noise, vt, to bs is given by:

The transfer functions of motor noise, vm, to tc and to t are given by:

The transfer functions of proprioceptive signal noise, vp, to tc and to t are given by:

The transfer functions of gravioceptive signal noise, vg, to tc and to t are given by:

The transfer functions of torque signal noise, vt, to tc and to t are given by:

Stimulus-dependent noise calculations required the calculation of the PSDs of torque, proprioceptive, and graviceptive signals evoked by the ss stimulus. For a linear system, the PSD of an output signal is given by the PSD of the input signal (ss in this case) multiplied by the squared absolute transfer function between the input and output (Bendat and Piersol 2000). The stimulus-evoked corrective torque PSD, Ptc, and stimulus-evoked torque sensor signal PSD, Pt, in the absence of noise are given by:

where Pss is the PSD of the surface-tilt stimulus. The PSD of the stimulus-evoked proprioceptive signal in the absence of noise is given by:

The PSD of the stimulus-evoked graviceptive signal, Pg, in the absence of noise is given by:

The mean-square value of a signal, E{()2}, can be calculated from the integral of the power spectral density, e.g. for the stimulus-evoked corrective torque:

where f is the increment between adjacent frequency points in the PSD.