## Abstract

A fundamental issue in motor control is how to determine the task goals for a given behavior. Here, we address this question by separately identifying the musculoskeletal and feedback components of the human postural control loop. Eighteen subjects were perturbed by two mechanical perturbations (gentle pulling from behind at waist and shoulder levels) and one sensory perturbation (movement of a virtual visual scene). Body kinematics was described by the leg and trunk segment angles in the sagittal plane. Muscle activations were described by ankle and hip EMG signals, with each EMG signal computed as a weighted sum of rectified EMG signals from multiple muscles at the given joint. The mechanical perturbations were used to identify feedback, defined as the mapping from the two segment angles to the two EMG signals. The sensory perturbation was used to estimate parameters in a mechanistic model of the plant, defined as the mapping from the two EMG signals to the two segment angles. Using the plant model and optimal control theory, we compared identified feedback to optimal feedback for a range of cost functions. Identified feedback was similar to feedback that stabilizes upright stance with near-minimum muscle activation, but was not consistent with feedback that substantially increases muscle activation to reduce movements of the body's center of mass or center of pressure. The results suggest that the common assumption of reducing sway may not apply to musculoskeletal systems that are inherently unstable.

## Introduction

A long-standing issue in human motor control is to determine the task goals for a given behavior and how those task goals are achieved by the nervous system. For example, reduction of center of mass (COM) or center of pressure (COP) movement is often assumed as the task goal for standing posture given the roles these variables play in postural stability. Identifying task goals is closely related to the degrees of freedom problem (Bernstein, 1967), since in a system with multiple mechanical degrees of freedom there may be multiple ways to achieve the same task goal (i.e., motor equivalence).

One approach to inferring task goals is to understand their relationship to feedback using optimal control theory. Here, we define feedback as the “neural control” component of the postural control loop (see Fig. 1). The theory is based on the assumption that feedback minimizes a cost function that penalizes control signals (muscle activations in this study) and state variables such as the COM or COP (Kuo, 1995, 2005; van der Kooij et al., 1999). The cost function serves to quantify the task goals for the motor behavior (Todorov, 2004). For example, if the task goal for standing were to reduce movements in the COM, then the cost function would penalize deviations of the COM away from its preferred value. The size of the COM penalty relative to the penalties on muscle activations would quantify the extent to which neural control is designed to reduce COM deviations at the expense of increasing muscle activation.

The power of the optimal-control approach is that it describes the relationship among the plant, feedback (see Fig. 1), and the cost function that describes the task goals. Given a model of the plant, for every cost function there is an optimal feedback that minimizes the cost function. Thus, by identifying the plant and feedback for a motor task such as standing, it is possible, in principle, to infer the cost function. The challenge is to identify the plant and feedback components of the closed-loop postural control system without opening the loop, which is not possible because the plant is unstable; if all sensory information is removed, the person falls. In this study, we used the joint input–output method (JIOM) of closed-loop system identification (Katayama, 2005). In the JIOM, one identifies the plant and feedback by measuring the responses of plant inputs (EMG signals) and plant outputs (body segment angles) to sensory and mechanical perturbations.

Fitzpatrick et al. (1996) were the first to apply the JIOM to postural control, using a single-joint (ankle) approximation of the body. Recently, we extended the method to identify the plant for a multijoint body (Kiemel et al., 2008). In the present study, we identified both the plant and feedback for a multijoint body, allowing us to infer the cost function. Our results show that neural feedback appears designed to nearly minimize the muscle activation necessary to stabilize upright stance. It does not produce the additional muscle activation that would be necessary to reduce COM or COP displacements.

## Materials and Methods

#### Experimental Methods

##### Subjects.

Eighteen subjects (nine males; nine females), aged 18–27, were recruited at the University of Maryland. The procedures used in the experiment were approved by the Institutional Review Board at the University of Maryland. All subjects received instructions for the test procedures. Informed written consent was obtained from all participants in the study. All the subjects were physically active, with no known musculoskeletal injuries or neurological disorders that might affect their ability to maintain balance.

##### Experiment setup.

Two weak continuous mechanical perturbations were applied to the subject from behind. A waist-level perturbation was applied using a spring with one end attached to a waist belt worn by the subject and the other end attached to a linear motor (LX80L; Parker Hannifin Corporation) placed directly behind the subject. Similarly, a shoulder-level perturbation was applied using a spring attached at one end to a shoulder strap worn by the subject and the other end attached to a different linear motor. The actual displacements of the motors were used as the mechanical perturbation signals. The waist and shoulder springs were weak with spring constants of 0.04 and 0.0157 N/mm, respectively. We used weak springs to reduce their stabilizing effect on the subject and thus reduce their effect on the plant.

A weak sensory perturbation was applied to the subject using a moving virtual visual scene. The subject stood surrounded by three screens (width, 3.05 m; height, 2.44 m; Fakespace), one in front of the subject and one on either side. Participants stood 107 cm from the front screen and centered between the two side screens. Visual displays were rear-projected to the screens at a frame rate of 60 Hz by JVC projectors (model DLA-M15U; Victor Company of Japan). CaveLib software (Fakespace) was used to generate a virtual moving visual scene consisting of three walls attached at right angles that coincided with the screens when the visual scene was not moving. Each wall consisted of 500 non-overlapping white small triangles (3.4 × 3.4 × 3.0 cm) with random positions and orientations on a black background. To reduce aliasing effects in the fovea region, no triangles were displayed on the front wall within a 30-cm-radius circular region directly in front of the participant's eyes. The display on each screen was varied with time to simulate rotation of the visual scene about the axis through the subject's ankles, assuming a fixed perspective point at the average position of the participant's eyes.

Subjects assumed a foot position with heels at a distance of ∼11% of their heights and feet pointing outward with an angle of 14° between each foot and the midline (McIlroy and Maki, 1997). The instruction to the subjects was to look straight ahead at the front screen and not to consciously resist any force from the waist belt and shoulder strap. There was only one experimental condition with all three perturbations applied simultaneously. There were 12 250 s trials for each subject (one subject had only 10 trials due to technical difficulties).

##### Kinematics.

The kinematics was captured by an Optotrak infrared position tracking system (Northern Digital). The shoulder (the scapula), hip (the greater trochanter), knee (the lateral femoral condyle), and ankle (the lateral malleolus) were measured by attaching four LED markers on the right side of the subject. The markers were sampled at 120 Hz. Only data from the ankle, hip, and shoulder markers were used in the present study.

##### EMG.

Muscle activity was measured using a multichannel telemetric surface EMG system (Zerowire; Aurion). Silver/silver chloride electrodes were placed on the right side of the body measuring activity from 12 muscles. Based on the two-joint approximation of the body described below, in the present study we only analyzed data from four muscles acting at the ankle (lateral gastrocnemius, medial gastrocnemius, soleus, and tibialis anterior) and three muscles acting at the hip (biceps femoris, rectus femoris, and semitendinosus). (Some of these muscles also act at the knee.) EMG signals were bandpass filtered between 10 and 1000 Hz and sampled at 2160 Hz.

##### Mechanical and visual perturbation signals.

We used filtered white-noise signals for perturbation signals, so that all three simultaneous perturbations would be statistically independent. White noise with a power spectral density of *P*_{0} was passed through a first-order filter with cutoff frequency of *f*_{c1} and an eighth-order Butterworth filter with a cutoff frequency of *f*_{c2}. The seed used to generate the white noise was different for every subject, trial, and perturbation signal. For the waist motor displacement, *P*_{0} = 4 cm^{2}/Hz, *f*_{c1} = 0.1 Hz, and *f*_{c2} = 5 Hz. For the shoulder motor displacement, *P*_{0} = 2.5 cm^{2}/Hz, *f*_{c1} = 0.1 Hz, and *f*_{c2} = 5 Hz. For visual-scene angle, *P*_{0} = 4.05 deg^{2}/Hz, *f*_{c1} = 0.02 Hz, and *f*_{c2} = 5 Hz. The initial and final 5 s of each 250 s signal were multiplied by increasing and decreasing ramps, respectively, to insure that the value of the signal at the beginning and end of the trial would be 0. Only the middle 240 s of each trial was analyzed.

#### Signal processing

##### Kinematic and EMG response signals.

The leg angle θ_{1}(*t*) and trunk angle θ_{2}(*t*) in the sagittal plane with respect to vertical were determined by the anterior–posterior and vertical displacement of the shoulder, hip, and ankle markers. Positive angles indicated forward lean. For the EMG activity of each muscle, the mean was subtracted from the raw EMG and then full-wave rectified, resulting in ankle EMG signals *u*_{11}(*t*), *u*_{12}(*t*), *u*_{13}(*t*), and *u*_{14}(*t*) and hip EMG signals *u*_{21}(*t*), *u*_{22}(*t*), and *u*_{23}(*t*).

##### Spectral analysis.

For any two signals *x*(*t*) and *y*(*t*), the power spectral densities (PSDs) *p _{xx}*(

*f*) and

*p*(

_{yy}*f*) and cross-spectral density

*p*(

_{xy}*f*), where

*f*is frequency, were computed using Welch's method (Bendat and Piersol, 2000) with 40 s Hanning windows and 50% overlap and then averaged across trials. The closed-loop frequency response function (FRF) from

*x*(

*t*) to

*y*(

*t*) is

*H*(

_{xy}*f*) =

*p*(

_{xy}*f*)/

*p*(

_{xx}*f*). Gain is the absolute value of

*H*(

_{xy}*f*) and phase is the argument of

*H*(

_{xy}*f*), converted to degrees. A positive phase indicates that

*y*(

*t*) was phase advanced relative to

*x*(

*t*).

##### Weighted EMG signals.

The four rectified ankle EMG signals *u*_{11}(*t*), *u*_{12}(*t*), *u*_{13}(*t*), and *u*_{14}(*t*) were used to compute the “weighted ankle EMG signal” *u*_{1}(*t*) = *w*_{11}*u*_{11}(*t*) + *w*_{12}*u*_{12}(*t*) + *w*_{13}*u*_{13}(*t*) + *w*_{14}*u*_{14}(*t*). Weights were adjusted using the MATLAB optimization toolbox to maximize the average coherence between the three perturbation signals and *u*_{1}(*t*) for frequencies from 0 to 5 Hz subject to the constraints that *w*_{1j} ≥ 0 for posterior muscles, *w*_{1j} ≤ 0 for anterior muscles, and |*w*_{11}| + |*w*_{12}| + |*w*_{13}| + |*w*_{14}| = 1. Similarly, the three hip EMG signals were used to define the “weighted hip EMG signal” *u*_{2}(*t*) = *w*_{21}*u*_{21}(*t*) + *w*_{22}*u*_{22}(*t*) + *w*_{23}*u*_{23}(*t*). Each weighted EMG signal was normalized by dividing by the square root of its power from 0 to 5 Hz, calculated by integrating its PSD.

##### Frequency binning.

Our use of a 40 s spectral window resulted in 200 frequency values from 0.025 to 5 Hz. To improve the accuracy of our FRF estimates, these 200 values were divided into 10 frequency bins, and PSD and cross-spectral density values were averaged within bins before computing FRFs. The lower frequency values of the bins were 0.025, 0.100, 0.200, 0.300, 0.425, 0.650, 0.975, 1.475, 2.225, and 3.325 Hz. These values were chosen to be approximately equally spaced on a log scale with additional adjustments at low frequencies to reduce large errors in FRF estimates. When plotting gains and phases as a function of frequency, we used the average frequency in each bin.

#### Identification of the feedback FRF

We used the joint input–output method of closed-loop system identification (Katayama, 2005; van der Kooij et al., 2005) to identify the feedback FRF based on the theoretical framework of Figure 1 describing postural control in the sagittal plane. We assumed a linear approximation for each process in the postural control feedback loop. Let *u*(*t*) be the 2 × 1 vector of muscle activation signals (one each for the ankle and hip), *y*(*t*) be the 2 × 1 vector of body segment angles (legs and trunk), *v*(*t*) be visual-scene angle, and *d*(*t*) be the 2 × 1 vector of motor positions, and let *U*(*f*), *Y*(*f*), *V*(*f*), and *D*(*f*) be their Fourier transforms, where *f* is frequency. The theoretical framework of Figure 1 can be described in the frequency domain by the following:
where *P*(*f*) is the open-loop FRF of the plant, *M*(*f*) is a FRF characterizing the effect of the mechanical perturbations on body segment angles, *N _{y}*(

*f*) is the Fourier transform of intrinsic noise in the plant,

*F*(

*f*) is the open-loop FRF of the feedback,

*S*(

*f*) is a FRF characterizing the effect the visual perturbation on the EMG signals,

*N*(

_{u}*f*) is the Fourier transform of intrinsic noise in the feedback. From Equation 1, it follows that closed-loop FRFs

*H*(

_{vy}*f*) and

*H*(

_{vu}*f*) describing kinematic and EMG responses to the visual-scene angle are related by the following: where we have used that the visual-scene angle is uncorrelated with the mechanical perturbations and intrinsic noise. Similarly, from Equation 2, it follows that the closed-loop FRFs

*H*(

_{du}*f*) and

*H*(

_{dy}*f*) describing EMG and kinematic responses to the mechanical perturbations are related by the following: Based on Equation 4, the feedback FRF was identified as the following: Here, it was critical that the number of mechanical perturbations equaled the number of segment angles so that the matrix

*H*(

_{dy}*f*) was invertible, assuming that the FRFs describing the responses to the two mechanical perturbations [the columns of

*H*(

_{dy}*f*)] were linearly independent.

When averaging across subjects, we chose to first average the closed-loop FRFs across subjects and then compute feedback as follows: Compared with identifying the plant and feedback for each subject and then averaging across subjects, this method reduced errors caused by subjects with low coherence between mechanical perturbations and response variables. For closed-loop FRFs and the feedback FRF, bootstrap SEs for log-gain and phase were computed using 1000 bootstrap resamples of the subjects (Zoubir and Boashash, 1998).

#### Fitted mechanistic plant model

We used empirical closed-loop responses to the visual perturbation to fit parameters in a mechanistic two-joint (ankle and hip) model of the plant. The model provides a mapping from the ankle and hip EMG signals *u*_{1}(*t*) and *u*_{2}(*t*) to the leg and trunk segment angles *y*_{1}(*t*) and *y*_{2}(*t*). The model is the same as that in the study by Kiemel et al. (2008), except that the previous model assumed a fixed ratio between *u*_{1}(*t*) and *u*_{2}(*t*), whereas in the present study we did not.

To summarize, the body was approximated by a two-joint inverted pendulum in the sagittal plane with dynamics linearized about vertical as follows:
The input variables were the ankle torque *T*_{1}(*t*) and the hip torque *T*_{2}(*t*), where anterior muscles produce positive torque. The output variables were the leg segment angle *y*_{1}(*t*) and the head, arms, and trunk (HAT) segment angle *y*_{2}(*t*), where a positive angle indicates a forward lean with respect to vertical. The parameter *g* = 9.81 m^{2}/s is the acceleration due to gravity. Parameters for body segment *k* (*k* = 1 for legs; *k* = 2 for HAT) are *m _{k}*, mass;

*l*, length;

_{k}*h*, height of the center of mass above the lower end of the segment; and

_{k}*J*, the moment of inertia about the lower end of the segment. Values for these parameters were computed based on the average mass

_{k}*m*= 70.8 kg and segment lengths of the subjects and anthropometric parameters from Winter (1990) as described by Kiemel et al. (2008):

*m*

_{1}= 20.7 kg,

*m*

_{2}= 48.0 kg,

*l*

_{1}= 0.826 m,

*h*

_{1}= 0.512 m,

*h*

_{2}= 0.301 m,

*J*

_{1}= 6.57 kg m

^{2}, and

*J*

_{2}= 7.07 kg m

^{2}.

The musculotendon actuator for joint *j* (*j* = 1 for ankle, *j* = 2 for hip) was modeled as follows:
where ϕ_{1}(*t*) = *y*_{1}(*t*) was the ankle angle, ϕ_{2}(*t*) = *y*_{2}(*t*) − *y*_{1}(*t*) was the hip angle, and *u*_{1}(*t*) and *u*_{2}(*t*) were the weighted ankle and hip EMG signals, respectively. Equation 8, which describes the EMG-to-torque mapping from *u _{j}*(

*t*) to the active joint torque

*T*

_{0j}(

*t*), was modeled as a second-order low-pass filter with angular eigenfrequency ω

_{0j}, damping fraction η

*, and DC gain γ*

_{j}*. The negative sign in the right-hand side reflects our sign convention that*

_{j}*u*(

*t*) > 0 corresponds to activation of posterior muscles. Total joint torque, the sum of active torque and passive torque, is produced by intrinsic stiffness

*k*and intrinsic damping α

_{j}*(Eq. 9).*

_{j}The 10 parameters ω_{01}, ω_{02}, η_{1}, η_{2}, γ_{1}, γ_{2}, *k*_{1}, *k*_{2}, α_{1}, and α_{2} of the plant model were estimated based on the mean empirical closed-loop FRFs *H̅ _{vu}(f)* and

*H̅*describing EMG and kinematic responses to the visual perturbation. Based on Equation 3, we computed the model parameters that minimized the errors

_{vy}(f)*H̅*−

_{vy}(f)*P(f)H̅*, where

_{vu}(f)*P*(

*f*) is the FRF of the plant model. Specifically, the MATLAB optimization toolbox was used to adjust parameters to minimize the objective function as follows: where

*f*is the average frequency for the

_{k}*k*th frequency bin,

*e*is the

_{j}*j*th standard basis vector, and “T” indicates transpose. All model parameters were constrained to be non-negative. We tried multiple sets of initial parameters to increase the chance of finding the global minimum.

#### Comparison of identified and optimal feedback

We used the fitted plant model to compute optimal feedback for various cost functions. First, we converted the plant model to state-space form **ẋ**(*t*) = **Ex**(*t*) + **Gu**(*t*), where **x**(*t*) = [y_{1}(*t*),*ẏ*_{1}(*t*), *y*_{2}(*t*),*ẏ*_{2}(*t*),*T*_{01}(*t*),*Ṫ*_{01}(*t*), *T*_{02}(*t*),*Ṫ*_{02}(*t*)]** ^{T}** is the vector of state variables,

**u**(

*t*) = [

*u*

_{1}(

*t*),

*u*

_{2}(

*t*)]

^{T}is the vector of muscle activation (EMG) signals,

**E**is an 8 × 8 matrix, and

**G**is a 8 × 2 column vector. To include a feedback time delay τ, we let

**u**(

*t*) =

**ũ**(

*t*− τ) and defined the modified plant from

**ũ**(

*t*) to

**x**(

*t*) as

**ẋ**(

*t*) =

**Ex**(

*t*) +

**Gũ**(

*t*− τ). To approximate the modified plant with a finite number of state variables, we converted the model to discrete time with a time step of δ = τ/15, leading to the discrete-time state-space model

**x**

_{k+1}=

**E**

_{τ}

**x**

_{k}+

**G**

_{τ}

**ũ**

_{k}, where

**x**

*=*

_{k}**x**(

*k*δ) with 30 additional state variables appended to the end and

**ũ**

*=*

_{k}**ũ**(

*k*δ).

Next, we defined the cost function Φ =½*E*[**x**_{k}^{T}**Ax**_{k} + **ũ**_{k}^{T}**Bũ*** _{k}*], where

**A**is a 38 × 38 positive semidefinite matrix and

**B**is a 2 × 2 diagonal matrix with

*B*

_{11}> 0 and

*B*

_{22}> 0, and “

*E*” stands for expected value. To compute the optimal feedback, we made two simplifying assumptions: (1) noise in the system is white and additive; and (2) sensory information based on the segment angles

*y*

_{1}(

*t*) and

*y*

_{2}(

*t*) provides accurate estimates of the state vector

**x**(

*t*). Using optimal control theory for a linear-quadratic regulator (Bryson and Ho, 1975), we computed the optimal control

**ũ**

_{k}=

**−Cx**

*that minimizes Φ, where*

_{k}**C**is a 2 × 38 matrix. To compute the optimal feedback FRF

**(**

*F̃***f*) from

**y**

*= [*

_{k}*x*,

_{k}_{1}(

*t*),

*x*,

_{k}_{3}(

*t*)]

^{T}to

**ũ**

_{k}, we considered a disturbance

**d**

*added to the right-hand side of the discrete-time version of equation (8) and then used equation (5) to compute*

_{k}*F̃**(

*f*). Finally, the optimal feedback FRF

*̃F**(

*f*) from

**y**(

*t*) to

**u**(

*t*) was computed as

*F**(

*f*) =

*e*

^{−2πifτ}

*F̃**(

*f*).

The above procedure computes the optimal feedback *F**(*f*) as a function of the cost-function matrices **A** and **B** and the feedback time delay τ. We considered a family of matrices **A** of the form **A** = **A**_{kin} + μ_{COP}^{2}**A**_{COP}. **A**_{kin} describes the penalties on displacements of kinematic variables, the leg and trunk angles and their velocities; μ_{COP}^{2}**A**_{COP} is the penalty on displacements of the COP. **A**_{kin} is zero except for its 4 × 4 upper-left submatrix, which we parameterized using the Cholesky decomposition **LL**^{T}, where **L** is a lower-triangular 4 × 4 matrix with 10 nonzero entries. A special case of **A**_{kin} is a penalty on the COM. For the plant model, the center of mass displacement is *b*_{1}*y*_{1} + *b*_{2}*y*_{2} where *b*_{1} = (*m*_{1}*h*_{1} + *m*_{2}*l*_{1})/*m* and *b*_{2} = *m*_{2}*h*_{2}/*m*. Therefore, a penalty of size μ_{COM} on the COM corresponds to an **L** with zero entries except for *L*_{11} = μ_{COM} *b*_{1} and *L*_{31} = μ_{COM} *b*_{2}. To define **A**_{COP}, we ignored the height of the ankle and used the small-angle approximation that COP position equals –*T*_{1}/*mg*, where *T*_{1} is ankle torque (Barin, 1992). In the plant model, *T*_{1} = *T*_{01} − *k*_{1}*y*_{1} − α_{1}*ẏ*_{1}. Therefore, **A**_{COP} = **aa**^{T}, where **a** is a vector with zero entries except for *a*_{1} = −*k*_{1}/*mg*, *a*_{2} = −α_{1}/*mg*, and *a*_{5} = 1/*mg*.

We adjusted the 10 parameters of **L** describing penalties on kinematic deviations, the penalty μ_{COP} on COP deviations, and the feedback time delay τ to minimize the error between the optimal feedback *F**(*f*) and the identified feedback *F*(*f*). For the fitting procedure, we fixed the penalty **B** on muscle activations to avoid degenerate cases in which *B*_{11} or *B*_{22} approach 0. Such degenerate cases violate the central idea of the optimal control approach that there are trade-offs between reducing variation in behavioral variables, such as the COM and COP, and reducing muscle activation. We chose the values *B*_{11} = 1 and *B*_{22} = 0.66, which yielded reasonable agreement between *F**(*f*) and *F*(*f*) for the fitted **L**, μ_{COP}, and τ. We defined the error between *F**(*f*) and *F*(*f*) as follows:
where ρ(*z*) = (Re *z*, Im *z*)^{T} for a complex number *z*, and *S _{jkl}* is the bootstrap variance–covariance matrix of ρ(

*F**

*−*

_{jk}(f_{l})*F*) based on 1000 resamples. We used the optimization toolbox of MATLAB to find the values of our 12 parameters that minimized the error. We tried multiple sets of initial parameters, all of which resulted in the same set of fitted parameter values. For the fitted cost function, entries of

_{jk}(f_{l})**L**were zero except for its first column, which we denoted by the 4 × 1 vector

**v**. Such an

**L**corresponds to a penalty on deviations of the kinematic variable

*y*

_{min}=

*v*

_{1}

*y*

_{1}+

*v*

_{2}

*ẏ*

_{1}+

*v*

_{3}

*y*

_{2}+

*v*

_{4}

*ẏ*

_{2}, where

*y*

_{1}and

*y*

_{2}are the leg and trunk angles. In Results, we describe the fitted cost function and time delay by the fitted values of

**v**, μ

_{COP}, and τ and their bootstrap SEs based on 100 resamples.

To quantify the amounts of sway and muscle activations for a given optimal feedback law, we inserted noise into the closed-loop discrete-time system **x**_{k+1} = (**E _{τ} − G_{τ}̃C)x**

_{k}+

**w**

*, where*

_{k}**w**

*is vector white noise with variance–covariance matrix δ*

_{k}**Q**. The variance–covariance matrix Cov[

**x**

*] is the symmetric positive-definite solution to the following: which we computed through iteration. We added noise to the joint torques in the right-hand sides of Equation 7, corresponding to a*

_{k}**Q**with zero entries except for

*Q*

_{33},

*Q*

_{44}, and

*Q*

_{34}=

*Q*

_{43}. We adjusted these noise parameters to match the empirical variance–covariance matrix of the leg and trunk angles. We added noise to joint torques, because it was the simplest way we found to match the empirical data. To quantity the amount of deviations in a sway variable, we used its root-mean-square (RMS) value computed from Cov[

**x**

*]. We quantified the amount of muscle activation as*

_{k}**ũ**

_{k}] =

**C**Cov[

**x**

_{k}]

**C**

^{T}.

## Results

Figure 2 shows example time series for 10 s from a single trial. Note that the ankle and hip weighted EMG signals have both positive and negative values (Fig. 2*D*,*E*). Based on our sign convention, the EMG signals are positive when posterior muscles are primarily active and negative when anterior muscles are primarily active. Also, note that the slower changes in the EMG signals appear to be positively correlated with leg and trunk angles (Fig. 2*F*,*G*).

### Identification of feedback

Closed-loop FRFs from mechanical perturbations to EMG signals and segment angles were used to identify the open-loop FRF describing feedback (see Materials and Methods, Comparison of identified and optimal feedback). Figure 3, *A* and *B*, shows the gains and phases of the mean closed-loop FRF *H̅ _{du}(f)* from the waist and shoulder mechanical perturbations to the ankle and hip EMG signals.

*H̅*has 2 inputs and 2 outputs so

_{du}(f)*H̅*is a 2 × 2 matrix at each frequency

_{du}(f)*f*. Similarly, Figure 3,

*C*and

*D*, shows the gains and phases of the mean closed-loop 2 × 2 FRF

*H̅*from the waist and shoulder mechanical perturbations to the leg and trunk segment angles.

_{dy}(f)It is difficult to mechanistically interpret the closed-loop FRFs *H̅ _{du}(f)* and

*H̅*, because they reflect the closed-loop interaction between the plant and feedback components of the postural control feedback loop as well as the properties of the mechanical perturbations (Fig. 1). However, the relationship between the two closed-loop FRFs depends only on the feedback component. Specifically, the inferred open-loop FRF of feedback is

_{dy}(f)*F(f)*=

*H̅*

_{du}(f)H̅_{dy}(f)^{−1}, whose gain and phase are shown in Figure 3,

*E*and

*F*. Feedback is the mapping from the leg and trunk segment angles to the ankle and hip EMG signals. At low frequencies, the gain of feedback is approximately constant (Fig. 3

*E*) and phase is near 0 (Fig. 3

*F*), indicating that the activations of ankle and hip muscles are approximately proportional to deviations of segment angles from their mean positions. As frequency increases, gains increase and phases also initially increase. This pattern is consistent with feedback that also depends on segment-angle velocities, that is, proportional-derivative (PD) control (Johansson et al., 1988; Peterka, 2000, 2002). With increasing frequency, phases reach a maximum between 94 and 103° at the seventh frequency bin (1.2 Hz) and then decrease. (The only departures from this overall pattern occur for the hip EMG components of feedback in the last frequency bin.) The eventual decrease in phase is consistent with the presence of a time delay in the feedback pathway (see Results, Comparisons of identified feedback to optimal feedback). The fact that the maximum phase advance reaches 90° even with a feedback time delay indicates that the identified feedback is not consistent with PD control; feedback also depends on segment-angle accelerations and perhaps higher derivatives.

### Fitted mechanistic model of the plant

Closed-loop FRFs from the visual perturbation to EMG signals and segment angles were used to fit parameters in a mechanistic model of the plant (see Materials and Methods, Fitted mechanistic plant model). Figure 4, *A* and *B*, shows the gain and phase of the mean closed-loop FRF *H̅ _{vu}(f)* from the visual-scene angle to the ankle and hip EMG signals.

*H̅*has 1 input and 2 outputs, so

_{vu}(f)*H̅*is a 2 × 1 column vector at each frequency

_{vu}(f)*f*. Similarly, Figure 4,

*C*and

*D*, shows the gain and phase of the mean closed-loop 2 × 1 FRF

*H̅*from the visual-scene angle to the leg and trunk segment angles.

_{vy}(f)Just as the relationship between closed-loop responses to mechanical perturbations depends on feedback, the relationship between the closed-responses to sensory perturbations depends on the plant. Specifically, *H̅ _{vy}(f)* =

*P(f)H̅*, where the plant

_{vu}(f)*P*(

*f*) is the 2 × 2 FRF describing the mapping from the ankle and hip EMG signals to the leg and trunk segment angles. Since we had a single sensory perturbation, we could not directly identify the plant. Instead, we computed

*P*(

*f*) by fitting parameters in a mechanistic plant model to minimize the errors

*H̅*−

_{vy}(f)*P(f)H̅*. The dashed lines in Figure 4,

_{vu}(f)*C*and

*D*, are the gains and phases of

*P(f)H̅*, which are close to the gains and phases of

_{vu}(f)*H̅*, indicating that the model fit is reasonably accurate. Parameter values of the fitted plant are given in the legend of Figure 4. A stability analysis (Kiemel et al., 2008) shows that the plant is unstable, as indicated by one unstable eigenvalue. Figure 4,

_{vy}(f)*E*and

*F*, shows the gains and phases of the fitted plant model

*P*(

*f*). Gains eventually decrease with increasing frequency, indicating that the plant model acts as a low-pass filter in response to muscle activation. Phase for the plant model approaches −180 at low frequencies for the hip EMG to trunk component of the plant. Phase for the other three components approach 0 at low frequencies. Note that these phase responses only describe the plant. Since the plant is unstable, it is not possible to directly observe these phase relationships in experiments.

### Comparisons of identified feedback to optimal feedback

We now address the question of whether the identified feedback based on experimental data (Fig. 3*E*,*F*) is similar to optimal feedback that minimizes some cost function (see Materials and Methods, Comparison of identified and optimal feedback). We considered a family of cost functions that penalizes kinematic deviations, COP deviations, and muscle activations. The kinematic penalties allowed for penalties on any combination of the leg and trunk angles and their velocities, which includes the COM as a special case. We assumed that feedback is optimized to minimize the cost function given a feedback time delay τ. We fixed the penalties on muscle activations and adjusted the kinematic penalties, the COP penalty, and τ to achieve the best fit of the optimal feedback with the identified feedback. The best fit was obtained with a feedback time delay τ of 127.8 ± 1.4 ms (fitted value ± bootstrap SE), no penalty on COP deviations (μ_{COP} = 0.000 ± 0.000 cm^{−1}), and a kinematic-penalty matrix of the form **A**_{kin} = **vv**^{T}. The kinematic penalty corresponds to a penalty on deviations of the kinematic variable *y*_{min} = *v*_{1}*y*_{1} + *v*_{2}*ẏ*_{1} + *v*_{3}*y*_{2} + *v*_{4}*ẏ*_{2}, where *y*_{1} and *y*_{2} are the leg and trunk angles and *v*_{1} = 20.7 ± 10.5 deg^{−1}, *v*_{2} = 10.5 ± 10.1 (deg/s)^{−1}, *v*_{3} = −38.7 ± 4.0 deg^{−1}, and *v*_{4} = −15.0 ± 14.4 (deg/s)^{−1}. Note that minimizing *y*_{min} is fundamentally different from minimizing the COM. Unlike the COM, which is minimized when the leg and trunk angles have opposite signs, *y*_{min} is minimized when the leg and trunk angles have the same sign. Thus, our fitted cost function is not consistent with feedback designed to reduce deviations in either the COM or COP. We will address the functional significance of the penalty of the fitted cost function on *y*_{min} below.

The gains and phases of the optimal feedback that minimizes the fitted cost function (Fig. 5*A*,*B*) are similar to identified feedback (Fig. 3*E*,*F*). Optimal feedback gains (Fig. 5*A*) have plateaus at low frequencies with values of the same order of magnitude as the corresponding plateaus of identified feedback gains (Fig. 3*E*). The increase of gain with increasing frequency is qualitatively similar to that of identified feedback, although the gains of optimal feedback increase at a higher rate. At high frequencies, gain was highest for the legs-to-ankle EMG component of feedback, followed by the legs-to-hip EMG component, the trunk-to-ankle EMG component, and finally the trunk-to-hip EMG component. This is the same order of gain values for identified feedback at higher frequencies. At lower frequencies, the order of gain values for identified feedback was less consistent, but gain was always highest for the legs-to-ankle EMG component, consistent with the optimal feedback. The phase of optimal feedback (Fig. 5*B*) shows nearly in-phase activation of ankle and hip muscles in response to rotation of either the leg or trunk segment, similar to identified feedback (Fig. 3*F*). Also, the dependence of phase on frequency is similar to that of identified feedback. Phase initially increases with increasing frequency, reaches a maximum slightly above 1 Hz, and then decreases due to the time delay.

To help understand the fitted optimal feedback, we compared it with optimal feedback for two special cases. Figure 5, *C* and *D*, shows that optimal feedback that minimizes muscle activation subject to the constraint that feedback must stabilize the plant. It was obtained by multiplying **A**_{kin} for the fitted cost function by 10^{−8}. Essentially, the same optimal feedback would be obtained with a small penalty on any sway variable; the effect of the small penalty is to insure that feedback stabilizes upright stance. If the penalties on all sway variables are zero, then the optimal feedback is zero and the system is not stable. Note that the minimum-activation feedback, like fitted feedback (Fig. 5*A*,*B*), is similar to identified feedback (Fig. 3*E*,*F*). Differences between the two examples of optimal feedback are primarily in the relative gains of the different feedback components, with fitted feedback more closely matching identified feedback. In contrast, optimal feedback designed primarily to reduce COM deviations (Fig. 5*E*,*F*) exhibits more dramatic differences to identified feedback. Its gains are an order of magnitude larger than those of identified feedback. Its hip phases are approximately antiphase to its ankle phases, unlike those of identified feedback.

Given that the FRFs of fitted optimal feedback and minimum-activation feedback are qualitatively similar in many respects, the question arises whether they are functionally similar. To address this question, we put noise into the closed-loop optimal feedback model in such a way to match the empirical variances and covariance of the leg and trunk angles (see Materials and Methods, Comparison of identified and optimal feedback). We then compared the amount of muscle activation required by the fitted optimal feedback and the minimum-activation feedback. Averaged across subjects, the variances of the leg and trunk angles were 0.360 and 0.525 deg^{2}, respectively, and their covariance was 0.168 deg^{2}. The closed-loop model with fitted optimal feedback reproduced this variance structure with noise parameters *Q*_{33} = 0.0205 N^{2}m^{2}deg^{2}s^{−1}, *Q*_{44} = 0.4315 N^{2}m^{2}deg^{2}s^{−1}, and *Q*_{34} = −0.0757 N^{2}m^{2}deg^{2}s^{−1}.

We fixed the noise levels in the model and computed the RMS of *y*_{min} and the amount of muscle activation as we varied the penalty on *y*_{min}. The result is illustrated by the continuous curve in Figure 6. The largest RMS of *y*_{min} and the smallest level of muscle activation occurs as the penalty approaches 0 (open circle), corresponding to the minimum-activation feedback. As the penalty on *y*_{min} increases, the RMS of *y*_{min} decreases as muscle activation increases. However, muscle activation initially increases slowly, so that the fitted feedback produces only 1.6 ± 0.6% more muscle activation than that necessary to stabilize upright stance while reducing the RMS of *y*_{min} by a factor of 1.9 ± 0.5 (filled circle). In contrast to *y*_{min}, a substantial reduction in either COM of COP variation requires a substantial increase in muscle activation (Fig. 6, dashed and dotted lines). These modeling results indicate that there is a range of feedback laws that stabilize upright stance with near-minimum muscle activation. The optimal feedback most similar to our empirically identified feedback lies in this range, whereas feedback that substantially reduces deviations in the COM or COP does not.

## Discussion

We used closed-loop system identification techniques to separately identify the properties of the plant and feedback within the postural control loop. By comparing identified feedback to optimal feedback for a range of cost functions, we found that identified feedback was similar to feedback that stabilizes upright stance while achieving near-minimum muscle activation. Identified feedback was not consistent with feedback that substantially increases muscle activation to reduce deviations in the COM or COP.

Note that our results indicate that feedback achieves only approximate minimization of muscle activation, not perfect minimization. We found that there is a range of feedback laws that achieve such approximate minimization, representing a type of motor equivalence (Bernstein, 1967). Out of the wide range of control laws that stabilize upright stance, there is a smaller but still substantial range of control laws that do so with near-minimum muscle activation. Under the experimental conditions of this study, the nervous system adopts a control law in this smaller range, an indication that there are no additional task goals that require substantially increasing muscle activation.

### Control of the COM and COP

It may seem surprising that the cost function most consistent with our data does not penalize displacements of the COM and COP, since control of these variables is necessary to stabilize upright stance. In particular, the COP and vertical projection of the COM onto the support surface must remain within the subject's base of support (Kuo and Zajac, 1993; Ito et al., 2004). To understand this issue it is helpful to distinguish between linear stabilization that ignores the nonlinear constraint resulting from a finite base of support (FBOS) and the usual notion of stabilization that requires that the COM and COP remain within the base of support. If feedback that linearly stabilizes upright stance leads to COM and COP deviations much less than the size of the base of support, then such feedback (for all practical purposes) also achieves FBOS stability. Our results suggest that this is the case for feedback in the postural control system under our experimental conditions. Specifically, feedback that linearly stabilizes upright stance with near-minimum muscle activation results in FBOS stability. Thus, there is no obvious functional advantage to further minimizing COM or COP displacements that would justify increasing muscle activation.

The preceding discussion illustrates that the cost function in an optimal control framework may not reflect all the task goals of the motor behavior, especially if the plant is unstable. This suggests a similar limitation in using uncontrolled manifold (UCM) analysis (Scholz and Schöner, 1999; Hsu et al., 2007) and the related minimum intervention principle (Todorov and Jordan, 2002; Valero-Cuevas et al., 2009) to understand such motor behaviors. UCM analysis is based on the idea that feedback control will primarily act to correct deviations in variables relevant to the task goals, resulting in a variance structure that reflects those task goals. However, if linearly stabilizing upright stance results in the COM and COP remaining within the base of support without the need for additional muscle activation, there is no reason to expect that the variance structure of the behavior will reflect the importance of controlling the COM and COP.

### Linear continuous feedback

The optimal feedback under the theoretical framework of this study is linear, a type of continuous feedback control. Recently, Asai et al. (2009) argued that intermittent control has an advantage over continuous control in that it provides more robust stabilization of upright stance in the presence of a feedback time delay. This conclusion was based on comparing intermittent control to PD continuous control, that is, feedback that only depends on the instantaneous values of the segment angles and their velocities. However, as we have pointed out (see Results, Identification of feedback), our identified feedback is inconsistent with PD control. Instead, identified feedback also depends on the accelerations and perhaps higher derivatives of the segment angles. The same is true of optimal feedback for the various scenarios considered in Figure 5. Thus, both Asai et al. (2009) and the current study indicate that continuous PD control has disadvantages if there is a time delay in the feedback pathway. Whether these disadvantages are overcome by using intermittent control or a more general continuous control is an open question. One piece of evidence in favor of continuous control is a study by van der Kooij and de Vlugt (2007). In a direct test of continuous versus intermittent feedback using platform perturbations, they found that most of the postural responses could be explained by continuous feedback.

Linear feedback is the optimal feedback for the control theory framework used in this study (Bryson and Ho, 1975). We approximated the postural control system as a linear-quadratic regulator (LQR), that is, the plant is linear and the cost function is quadratic. The assumption that the plant is linear is a plausible approximation under the conditions of the current study: healthy individuals standing on a noncompliant fixed surface subjected to weak perturbations. Under these conditions, the displacements of the COM and COP are relatively small compared with the base of support provided by the feet.

Our use of LQR theory also assumes that the perturbations we used to identity feedback and the plant do not cause the nervous system to change the feedback law. The closed-loop responses to visual perturbations, which reflect properties of both feedback and the plant, in this study are similar to those in a previous study without mechanical perturbations (Kiemel et al., 2008), suggesting that mechanical perturbation do not substantially change feedback. The question of whether our visual perturbation changed feedback is more complex due to sensory reweighting (Peterka and Benolken, 1995; Oie et al., 2002; Peterka, 2002). Visual-scene motion is known to cause the nervous system to “down-weight” vision and compensate by “up-weighting” other sensory inputs, such as those from proprioception and the vestibular system. Such compensation suggests that feedback, which reflects integration by the nervous system of sensory inputs from all modalities, may not be strongly affected by our small visual perturbation. We are currently using the techniques of the present study to test this hypothesis.

Although linear approximations may be appropriate for the experimental conditions of the current study, this is not necessarily true under more challenging conditions that change the pattern of ankle versus hip EMG activity (Horak and Nashner, 1986; Horak et al., 1990; Park et al., 2004). These changes may be due to the nonlinearity of plant due to the finite base of support discussed above. One approach to approximating optimal control with a finite base of support is to change the cost function to more strongly penalize large deviations of the COM or COP (Li and Levine, 2009). Optimal feedback with such a cost function is no longer a linear function of segment angles and their derivatives.

### Single-joint versus multijoint models of posture

Using a multijoint model to identify the plant and feedback allows one to more easily distinguish among different hypotheses of neural control. For example, the phase difference between the ankle and hip EMG outputs of feedback was the clearest distinction between feedback designed to minimize muscle activation (Fig. 5*C*,*D*) and feedback designed to reduce COM displacements (Fig. 5*E*,*F*). Such a distinction would, of course, not exist for a single-joint model.

Also, a multijoint model can substantially improve the accuracy of identified feedback. For example, Fitzpatrick at al. (1996) used a single mechanical perturbation to identify feedback from a single body segment (shank) to a single EMG signal (soleus). The phase of their identified feedback monotonically increased with frequency and did not show the decrease at higher frequencies that would be expected given a feedback time delay (van der Kooij et al., 2005; Fig. 5*B*,*D*,*F*). One possibility is that part of the ankle muscle activation that Fitzpatrick et al. attributed to movement of the lower body was actually due to movement of the upper body. At higher frequencies, the movements of the lower and upper body can be substantially out of phase (Creath et al., 2005; Fig. 3*D*). By allowing ankle muscle activation to depend on movements of both the lower and upper body, we avoided this potential misattribution and identified feedback whose phase did show the expected phase decrease at higher frequencies (Fig. 3*F*).

### Conclusions

The feedback identified in this study suggests that, for unperturbed or weakly perturbed stance, the nervous system is primarily concerned with minimizing muscle activation. The nervous system does not produce substantially more activation than that necessary to stabilize upright stance. In particular, feedback does not seem designed to minimize COM or COP displacements. It will be important to determine whether this strategy remains when stability is challenged by conditions that increase postural sway or reduce the base of support.

## Footnotes

This work was supported by NIH Grant R01 NS35070 (J.J.J., principal investigator) and by a doctoral dissertation award from the Neuroscience and Cognitive Science Graduate Program at the University of Maryland (Y.Z.).

- Correspondence should be addressed to Tim Kiemel, Department of Kinesiology, University of Maryland, College Park, MD 20742-2611. kiemel{at}umd.edu