## Abstract

In the literature, it has been hotly debated whether the brain uses internal models or equilibrium point (EP) control to generate arm movements. EP control involves specification of EP trajectories, time series of arm configurations in which internal forces and external forces are in equilibrium; if the arm is not in a specified EP, it is driven toward this EP by muscle forces arising due to central drive, reflexes, and muscle mechanics. EP control has been refuted by researchers claiming that EP trajectories underlying movements of subjects were complex. These researchers used an approach that involves applying force perturbations during movements of subjects and fitting a mass-spring-damper model to the kinematic responses, and then reconstructing the EP trajectory using the estimated stiffness, damping, and measured kinematics. In this study, we examined the validity of this approach using an EP-controlled musculoskeletal model of the arm. We used the latter model to simulate unperturbed and perturbed maximally fast movements and optimized the parameter values of a mass-spring-damper model to make it reproduce as best as possible the kinematic responses. It was shown that estimated stiffness not only depended on the “true” stiffness of the musculoskeletal model but on all of its dynamical parameters. Furthermore it was shown that reconstructed EP trajectories were in agreement with those presented in the literature but did not resemble the simple EP trajectories that had been used to generate the movement of the model. It was concluded that the refutation of EP control on the basis of results obtained with mass-spring-damper models was unjust.

## INTRODUCTION

Theories proposed for the control of goal-directed (arm) movements come in two types: internal model (IM) control theories and equilibrium-point (EP) control theories. IM control theories rely on internal models of the dynamics of the musculoskeletal system to generate the muscle stimulation patterns (e.g., Kawato 1999; Mehta and Schaal 2002; Schweighofer et al. 1998; Shidara et al. 1993; Todorov and Jordan 2002; see for a comprehensive overview of the different types of IM controllers Wolpert et al. 1998). EP control involves the specification of an arm configuration in which internal forces and external forces are at equilibrium, or an EP trajectory, i.e., a time series of such configurations (e.g., Feldman et al. 1990; Gribble et al. 1998; McIntyre and Bizzi 1996). According to EP control, muscle forces are not explicitly computed but rather arise when the limb is not in the specified equilibrium configuration, due to central drive, reflexes and muscle mechanics. At least for single-joint movements, under EP control, there is no need for an internal dynamics model of the musculoskeletal system; only a mapping from the neural inputs to the muscles to the equilibrium arm configurations and stiffness is required (Kistemaker et al. 2006).

Although EP control is parsimonious and allows for a natural integration of the control of posture and the control of movement (Ostry and Feldman 2003), it has been rejected by many authors after Gomi and Kawato (1996, 1997) had estimated joint stiffness during fast arm movements of human subjects. Gomi and Kawato (1996, 1997, see also Katayama and Kawato 1993; Popescu et al. 2003) argued that under EP control, the net moments driving the arm are the product of the stiffness and the difference between the actual movement trajectory and the equilibrium-point trajectory. To reconstruct the EP trajectory, stiffness K (and damping B) were first estimated in these studies in vivo by subjecting the human controlled musculoskeletal system to perturbations. The parameter values *K* (stiffness), *B* (damping), and *I* (inertia) of the second-order mass-spring-damper model (that will be referred to as the KBI-model from here on) are optimized to achieve a best fit between the experimentally observed perturbation responses and the KBI-model's responses. Then estimated stiffness, damping, and inertia and the measured kinematics are used to calculate the EPs. Gomi and Kawato (1996, 1997) reconstructed EP trajectories using such a KBI approach and concluded that the EP trajectories were not “simple,” i.e., they did not resemble the actual movement trajectories. A typical reconstructed EP trajectory first led the actual trajectory to generate the accelerating moment and then fell behind the actual trajectory to generate the decelerating moment. Furthermore, it had a velocity profile with multiple peaks, which was very different from the actual velocity profile with only one distinct peak. Obviously the calculation of complicated EP trajectories by the CNS would require a model of the dynamics of the musculoskeletal system and would therefore obliterate the computational attractiveness of EP control (e.g., Gomi and Kawato 1996; Wolpert and Ghahramani 2000).

In a recent study (Kistemaker et al. 2006), we used a musculoskeletal model of the arm to explore the feasibility of EP control for fast arm movements. This model contains a substantial amount of biological detail and in particular contains the elements and characteristics that make life difficult for any control theory (admittedly, the model does not contain individual motor units, but given the size principle this does not make it fundamentally easier to control). With an EP controller that combined open- and closed-loop control and using simple constant-speed EP trajectories, we could make our musculoskeletal model reproduce very well experimentally observed maximally fast elbow flexion and extension movements. According to the arguments of Gomi and Kawato (1996, 1997), fast movements can only be generated using simple EP trajectories if the stiffness is high. However, the stiffness in our model, more precisely the intrinsic low-frequency elbow joint stiffness (*K*_{ilf}; defined here as the change in steady-state elbow joint moment per unit change in steady-state joint angle), fell low in the range of experimentally estimated stiffness values reported in the literature (Kistemaker et al. 2007). This calls for a revisiting of the estimation of stiffness and the reconstruction of equilibrium trajectories.

The purpose of the present study was to examine the validity of using KBI models for estimation of stiffness and damping of the musculoskeletal system and for reconstruction of EP trajectories. It needs no argument that this can only be done if the true stiffness and damping of the musculoskeletal system and the true EP trajectory are known. For that reason, we employed a modeling and simulation approach, using the EP-controlled musculoskeletal model that has recently been shown to produce realistic single-joint movements (see Kistemaker et al. 2006). We first used the model to simulate unperturbed and perturbed movements at three different levels of open-loop intrinsic stiffness (*K*_{ilf}). Subsequently, we optimized the parameter values of a KBI model to make it reproduce as best as possible the responses of our controlled musculoskeletal model to perturbations, and we compared the parameter values so obtained to the “true” parameters of the controlled musculoskeletal model. Second, we used the estimated stiffness and damping values to reconstruct the EP trajectory and compared this trajectory with the EP trajectory that had served as an input for the EP controller.

## METHODS

### Outline of the study

Simulations were carried out with a musculoskeletal model of the arm controlled by a hybrid open- and closed-loop EP-controller. In the present study, we used the musculoskeletal model and the EP controller with simple constant-speed EP trajectories to simulate elbow extension movements from 120 to 60° at three different levels of intrinsic low-frequency joint stiffness (i.e., *K*_{ilf} = 16, 10 and 5 Nm·rad^{−1}), obtained by manipulating the level of co-contraction/open loop neural input to the muscle (STIM_{open}). For each stiffness level, we simulated the responses to moment perturbations at three different onsets depending on the elbow joint position. In line with the approach advocated in the literature (e.g., Bennet et al. 1992; Gomi and Kawato 1996, 1997; Gomi and Osu 1998; Popescu et al. 2003), these responses were approximated by a KBI model; the parameter values of the KBI model were optimized to achieve a best fit between the responses of the controlled musculoskeletal model to the perturbations and the KBI model's responses to these same perturbations. Subsequently, the parameter values of the KBI model so obtained were used in an attempt to reconstruct the EP trajectory that had served as an input for the EP-controlled musculoskeletal model. In the following text, the musculoskeletal model, the EP controller and the simulation procedures are described in detail.

### Musculoskeletal model of the arm

The musculoskeletal model of the arm used in this study (Fig. 1) has been described in full detail elsewhere (Kistemaker et al. 2006; see also Van Soest and Bobbert 1993). It consisted of three rigid segments (representing forearm, upper arm, and shoulder blade), interconnected by two hinge joints (representing glenohumeral joint and elbow joint). Only elbow flexion/extension movements in the horizontal plane were allowed. The lower arm was actuated by four lumped muscles: a monoarticular elbow flexor (MEF; representing m. brachioradialis, m. brachialis, m. pronator teres, m. extensor carpi radialis), a monoarticular elbow extensor (MEE; representing m. triceps brachii caput laterale, m. triceps brachii caput mediale, m. anconeus, m. extensor carpi ulnaris), a biarticular elbow flexor (BEF; representing m. biceps brachii caput longum and caput breve) and a biarticular elbow extensor (BEE; m. triceps brachii caput longum). The muscles were modeled as Hill-type units consisting of a contractile element (CE), a parallel elastic element (PE), and a series elastic element (SE). Activation dynamics, describing the relation between neural input to the muscle and active state was modeled according to Hatze (1981; see also Kistemaker et al. 2005). In line with Hatze's description, the neural input to the muscle will be termed STIM throughout this study, and active state will be referred to as *q*. A variable step-size ODE solver based on the Runge-Kutta (4, 5) formula was used to numerically solve the differential equations of the musculoskeletal model.

### EP controller

The EP trajectory was defined as a “ramp” trajectory from the initial (120°) to the final (60°) position. From the EP trajectory, desired CE length (λ), and desired CE contraction velocity (λ̇) were derived, as outlined in the following text. The total neural input to the muscle (STIM_{h}) generated by the hybrid open- and closed-loop EP controller equaled with STIM_{h} is calculated by: (1)

The expression {*x*}_{0}^{1} means that values of *x* higher than 1 were set to 1 and negative values were set to 0.

The STIM_{open} is a vector that refers to the open-loop part of the neural inputs to the modeled muscles [STIM_{MEE_open}; STIM_{MEF_open}; STIM_{BEF_open}; STIM_{BEE_open}], i.e., the neural inputs to the muscles that, in the absence of external forces, bring the arm to the desired EP on the basis of intrinsic muscle properties even if no neural feedback is operative. In a previous study, it was shown that *1*) EPs could be set open-loop over the whole physiological range of motion, *2*) each EP could be set with numerous different STIM_{open}, with each yielding a particular *K*_{ilf} (Kistemaker et al. 2007). *K*_{ilf} was defined as the change in steady-state elbow joint moment per unit change in steady-state joint angle. λ is a vector with the desired steady state CE lengths of the modeled muscles [λ_{MEE}; λ_{MEF}; λ_{BEE}; λ_{BEF}] and λ̇(t) is the time derivative of λ(*t*). *k*_{p} and *k*_{d} denote the feedback constants (see following text). Feedback of CE length (*l*_{CE}) and CE contraction velocity (*v*_{CE}) was assumed to be linear and delayed by 25 ms; the time delay in the feedback loop (δ) was implemented using a fifth-order Padé approximation (Golub and Van Loan 1989). The procedure of calculating STIM_{open} and λ for a given EP at a desired level of *K*_{ilf} involved the following steps: *1*) calculating muscle-tendon complex lengths (*l*_{MTC}) in the EP using the relation between *l*_{MTC} and shoulder and elbow angle (as described in Kistemaker et al. 2006). *2*) Calculating isometric CE forces (*F*_{CE_isom}) for all muscles such that: *a*) the sum of the moments exerted by these muscles relative to the elbow joint axis equals zero (Σ*M*_{MUS} = 0); *b*) the sum of the low-frequency stiffnesses of the muscle-tendon complexes (under open loop control) at the EP considered resulted in the desired low-frequency elbow joint stiffness *K*_{ilf} (i.e., the partial derivative of *M*_{MUS} with respect to the elbow joint angle at the EP considered); and *c*) sum of *F*_{CE_isom} is as small as possible. And *3*) calculating the steady-state CE lengths at the EP considered from *F*_{CE_isom}.

*Criterion c* in *step 2* was incorporated because several STIM_{open} vectors exist that yield a certain level of *K*_{ilf} in the desired EP. To deal with this indeterminacy, the STIM_{open} was selected yielding the smallest amount of total isometric CE force which is likely to minimize metabolic demands. Steps 1–3 are repeated for the whole EP trajectory in time steps of 0.01 s for all three desired stiffness values, i.e., 16, 10, and 5 Nm·rad^{−1}. λ̇, the vector of the desired CE contraction velocities, was calculated by dividing the difference between two successive λs by 0.01 s. All in all this yields three different time histories of STIM_{open}(*t*), λ(*t*), and λ̇ (*t*).

To facilitate the comparison between true and reconstructed EP trajectory, continuous control was used throughout this study, even though in a previous study, it was shown that when open-loop components were sent out intermittently, maximal attainable movement speed increased (Kistemaker et al. 2006).

### Optimization of feedback gains and duration of the EP trajectory

For each *K*_{ilf} level, we used a grid search to find a combination of duration of the EP trajectory, *k*_{p} and *k*_{d} that led to a movement with the highest φ̇_{max} and with an RMS error of <2° in the elbow angle time history. This RMS error value was based on a previous study in which similar simulation results were compared with experimental data (Kistemaker et al. 2006).

### Perturbations and fitting procedure

Moment perturbations were applied at the instant that the elbow joint reached one of three different angles: early (100°), mid (80°), and end (60°). Following Popescu et al. (2003), we used “square-wave” moment perturbations: the perturbing moment instantaneously switched from 0 to 20 Nm to stay there for 15 ms, subsequently switched instantaneously to −20 Nm to stay there for 15 ms, and then switched back to zero. The responses of the musculoskeletal system (Φ_{A}(*t*)) to the perturbations were obtained by subtracting the perturbed movements (φ_{P}(*t*)) from the unperturbed movements (φ_{U}(*t*)) (2)

In line with the literature (e.g., Popescu et al. 2003), the values for stiffness (*K*) and damping (*B*) of a second-order linear KBI model were optimized such that during the first 0.05 s after perturbation onset, the sum of the squared difference between Φ_{A}(*t*) and simulated response of the KBI model (Φ_{KBI}(*t*)) was minimized. The value for inertia (*I*) was set to be identical to that of the musculoskeletal model (i.e., 0.078 kg · m^{2}). The optimal values for *K* and *B* were identified using a Nelder-Mead simplex search method (Lagarias et al. 1998). This approach will henceforth be referred to as the KBI-approach to estimating stiffness.

### Reconstruction of the EP trajectory

The EP trajectory (φ_{eq}(*t*)) was reconstructed according to the approach proposed by Gomi and Kawato (1996, 1997). First, the controlled musculoskeletal system was simplified to a second-order linear model (3) where *M*_{KBI} is the moment delivered by the KBI model and *K*(*t*), *B*(*t*), and *I* are the stiffness, damping and inertia parameters of the KBI model. *K*(*t*) and *B*(*t*) were obtained by linear interpolation of *K* and *B* estimated at positions early, mid, and end. From *Eq. 3*, the EP trajectory φ_{eq}(*t*) can be reconstructed (4)

It might be argued that in reconstructing φ_{eq}(*t*) one should not use the absolute velocity, but rather a reference velocity because in the control scheme of the EP-controller damping is also relative to a reference velocity. That is, instead of *Eq. 3* one should use (5)

For sufficiently small Δ*t* (0.001 s was used throughout this study), φ̇_{eq}(*t*) can adequately be described by (6) so that the expression for calculating φ_{eq}(*t*) becomes (7)

At *t* = 0, the system is in equilibrium, hence φ_{eq}(0) = φ_{U}(0) and φ̇_{eq} (0) = 0. This approach, either as described by *Eqs. 4* or *7*, will henceforth be referred to as the KBI approach to estimating EP trajectories.

## RESULTS

Figures 2 *A* and 3 show the perturbed and unperturbed simulated movements for all three stiffness conditions for all three perturbation instants (see also Table 1). The unperturbed movements produced by the musculoskeletal model were qualitatively similar to experimental data reported in a previous study (Kistemaker et al. 2006). In that study, it was shown that the EP controller was capable of making the musculoskeletal model reproduce experimentally observed kinematic data. Figure 4 shows the perturbation responses of the model and experimentally observed responses during similar single-joint arm movements by Popescu et al. (2003). The responses of the musculoskeletal model to the perturbations were somewhat faster than responses observed in human subjects. This is due to the fact that in our simulations the perturbing moment was a square wave, whereas in the study of Popescu et al. (2003), it was intended to be a square wave but in reality increased more gradually. Nonetheless, the responses of the musculoskeletal model to the perturbations resemble those of human subjects (Popescu et al. 2003; see Fig. 4). These results lend further support for our contention that the controlled musculoskeletal model captures the salient features of the real system.

Figure 2 (*B–D*) shows the elbow angle responses to the moment perturbations, obtained by subtracting the perturbed from the unperturbed movement as well as the simulated responses of the KBI model, for the 16 Nm · rad^{−1} condition. The simulated responses when *K*_{ilf} was set to 10 and 5 Nm · rad^{−1} were quantitatively similar (data not shown). As can be appreciated from Fig. 2, the responses of the KBI model with optimized parameter values were very similar to those of the musculoskeletal model. This finding is in accordance with previous studies in which it was found that the response of the musculoskeletal system to small perturbations can be adequately approximated with a second order KBI model (e.g., Agarwal and Gottlieb 1977; Hunter and Kearney 1982; Winters and Stark 1988). Table 2 presents the values of *K* and *B* estimated using the KBI approach for the different stiffness conditions as well as the true stiffness values of the musculoskeletal model (*K*_{ilf}). For all stiffness conditions and all perturbation onsets, the estimated stiffness values were substantially higher than the true stiffness values.

The difference between estimated and true stiffness values was greatest near the end of the movement; estimated stiffness was on average almost tree times as high as the true stiffness. This result can be understood from the series connection between an elastic SE and visco-elastic CE. The slope of the force-velocity relationship, i.e., the damping, becomes less with increasing speed. Near the end phase of the movement, *v*_{CE} was lower than at the early and mid phases, and hence CE damping was higher. With higher damping, the moment perturbation gave rise to a smaller change in CE length, and as the change in SE length depends only on the magnitude of the force step, the total muscle length change was smaller in the end phase. Consequently, the stiffness, being the change in muscle force (moment) divided by the change in length (joint angle), estimated using the KBI approach is higher in the end phase than in the other phases.

The estimated values for *K* and *B* were used to reconstruct EP trajectories according to the method proposed by Gomi and Kawato (1996, 1997) (see *Eq. 4*). We have also used a similar method that takes into account a reference velocity (see *Eq. 7*). The EP trajectories so obtained are shown in Fig. 5*A*, together with the true EP trajectory. In agreement with EP-trajectories reconstructed in experimental studies (e.g., Gomi and Kawato 1996, 1997), our EP trajectories reconstructed using the KBI approach were N-shaped, i.e., they first lagged and then led the true EP trajectory. However, the reconstructed EP trajectories showed poor resemblance with the true EP trajectory, especially during the second half of the movement. It made little difference whether or not a reference velocity was taken into account (Fig. 5*A*). One might argue that the stiffness values estimated would be higher when fitting the response over a longer period than 50 ms as reflexes from muscle spindles take part in the response. To examine the influence of higher estimated stiffness, we reconstructed the EP trajectory using stiffness values up to 10 times as high as those estimated (Fig. 5*B*). It was found, as suggested in the literature (e.g., Gomi and Kawato 1996, 1997; Katayama and Kawato 1993; Popescu et al. 2003; Wolpert et al. 1998) that the increased stiffness caused the reconstructed EP trajectories to increasingly better approximate the actual movement. However, the reconstructed EP trajectories never resembled the true EP trajectory that had been used to simulate the movements (see Fig. 5*B*).

## DISCUSSION

The purpose of the present study was to examine the validity of using KBI models for estimation of the stiffness of the musculoskeletal system (stiffness was defined in this study as the change in steady-state elbow joint moment per unit change in steady-state joint angle) and for reconstruction of EP trajectories. Because this required the true stiffness and EP trajectory to be known, we employed a modeling and simulation approach. We used an EP controlled musculoskeletal model to simulate unperturbed and perturbed fast elbow extension movements. We then optimized the stiffness and damping of a KBI model to approximate as best as possible the response of the controlled musculoskeletal model to the perturbations and found that the obtained stiffness did not match the true stiffness of the musculoskeletal model. Next, the KBI model with optimized stiffness and damping was used to reconstruct the EP trajectory, and it was found that this reconstructed trajectory did not resemble the true EP trajectory that had served as an input for the EP controller to generate the movement. In the following text, we shall elaborate on these findings.

The stiffness values estimated using the KBI approach were substantially higher than the true stiffness values calculated directly from the model. To make sure that the implementation of the estimation procedure was not responsible for the difference, we repeated our simulations of unperturbed and perturbed movements using an actual KBI model (see *Eq. 5*) instead of our musculoskeletal system; reassuringly, applying the estimation procedure yielded exactly the stiffness and damping of the KBI model used in the simulations. In theory, a poor fit of the KBI model could explain the difference between estimated stiffness and true stiffness of the musculoskeletal model. However, the responses of the KBI model were nearly identical to those of the musculoskeletal model, so a poor fit of the KBI model can be ruled out as a possible explanation. Another possible explanation for the discrepancies in stiffness values could be related to the way in which stiffness was calculated from the musculoskeletal model. The true stiffness was calculated directly from the uncontrolled musculoskeletal model (see Kistemaker et al. 2007), i.e., without the influence of feedback of muscle (CE) length and contraction velocity, whereas the parameters of the KBI model were fitted essentially to the total dynamics of the controlled musculoskeletal model. However, only the first 50 ms of the responses after the start of the perturbation were used in the fitting of the KBI model. Because several time delays exist in the model (i.e., feedback delay and delays due to activation and contraction dynamics; see Kistemaker et al. 2006), the contribution of feedback to the stiffness and damping during this 50-ms window may be assumed to be negligible. As a matter of fact, such short time windows have been used in the literature for the very purpose of minimizing the role of feedback from muscle spindles (e.g., Popescu et al. 2003). To verify the assumption, we simulated perturbations with and without feedback and found the responses to be nearly identical (results not presented); applying the estimation procedure (*Eq. 5*) yielded a maximal difference in estimated stiffness of 1 Nm · rad^{−1}. To conclude, the difference between estimated stiffness and true stiffness of the musculoskeletal model was not caused by flaws in the estimation procedure.

Let us try to explain why the stiffness estimated using a KBI-approach is different from the true stiffness of the musculoskeletal model. We suspect that the explanation lies in the neglect of higher-order components of a controlled musculoskeletal model when its behavior is approximated with a second-order KBI model. A KBI model is a second-order mechanical system, as the skeleton is actuated by damper and spring that in essence work in parallel. In the real system, however, the skeleton is actuated by muscles that have a visco-elastic CE in series with an elastic tendon (tendon is used here as a shorthand for all tendinous tissue “seen” by the CE). As a consequence, and in contrast with the KBI model, angular velocity and CE contraction velocity are not directly related. Such a muscle is at least of order one and in combination with a 1-DF mechanical system, this yields a third- or higher-order dynamical system. Due to the differences in order, the parameter values of the second-order KBI model used to approximate the perturbation response of the higher-order real system will depend on the frequency content of that perturbation. When low-frequency perturbations are applied to the musculoskeletal model, the joint stiffness estimated by fitting a KBI model resembles the stiffness of the controlled musculoskeletal model. For example, if a constant moment is applied, the stiffness estimated equals the change in moment divided by the observed change in steady-state position and resembles the stiffness of the controlled musculoskeletal model. Unfortunately, low-frequency perturbations cannot be used in human subject experiments because supraspinal systems will take part in the response to the perturbation. For that reason, experimenters use fast perturbations. However, when fast perturbations are applied to the musculoskeletal model, the viscous behavior of CE “prevents” rapid changes in CE length. As a result, the length changes of the muscle-tendon complex are mostly due to changes in tendon length. Consequently, joint stiffness estimated during fast perturbations is determined mostly by tendon stiffness, which is substantially higher than the true stiffness of the muscle-tendon complex. In fact, very fast perturbations, known as the “quick release method,” are used to estimate tendon stiffness in vivo (e.g., De Zee and Voigt 2001; Hof 1999). In sum, the stiffness estimated using KBI models is not only dependent on the stiffness of the musculoskeletal system but on all its dynamical parameters, the contribution of which depends on the nature of the perturbation. Differences in the nature of the perturbation might (partly) explain why different authors testing subjects in the same condition arrive at different stiffness estimates. For example, when testing stiffness in a “do-not-resist” condition, Laquaniti et al. (1982) arrived at elbow joint stiffness values ≤33 Nm · rad^{−1}, while Popescu et al. 2002 arrived at values of only 5 Nm · rad^{−1}.

In the preceding text, we have established that the stiffness and damping estimated using the KBI approach depends on the nature of the perturbation. Furthermore we have shown with our model that the stiffness estimated using a KBI approach will not represent the true stiffness of the musculoskeletal system. When it comes to reconstructing EP trajectories, one might argue that it does not matter whether the stiffness and damping estimated depends on the nature of the perturbation as long as they reflect the actual dynamical properties of the system for a given situation. However, our simulation results show that the EP trajectory reconstructed using the estimated stiffness and damping did in fact not resemble the true EP trajectory that had served as input for the EP-controlled musculoskeletal model. (Again, this was not an implementation problem: we repeated our simulations of unperturbed and perturbed movements using an actual KBI model instead of our musculoskeletal system and found the reconstructed EP trajectory to be identical to that which had served as input for the EP-controlled KBI model.) The reconstructed trajectory was N-shaped, whereas the true EP trajectory was not (Fig. 5*A*). At first glance, one would think that this is due to the estimated stiffness being lower than the true stiffness of the musculoskeletal model, but the opposite was true. In fact, reconstructing the EP trajectory using higher stiffness values only caused the reconstructed trajectory to better approximate the actual movement (Fig. 5*B*) in accordance with the generally accepted notion (e.g., Gomi and Kawato 1996, 1997; Katayama and Kawato 1993; Popescu et al. 2003; Wolpert et al. 1998). However, we were not interested in reconstructing the actual movement, our challenge was to reconstruct the true EP trajectory.

The main reason why the reconstructed EP trajectory did not match the true EP trajectory is that in the KBI approach, an EP is reconstructed by calculating the difference “needed” between the EP and actual position (and velocity) to make the KBI model generate a moment equal to the measured net moment (e.g., Gomi and Kawato 1996, 1997; Katayama and Kawato 1993;Popescu et al. 2003). However, like in most other suggested EP-control schemes, the difference between EP and actual position (and the difference between reference velocity and actual velocity) determines the *neural input* to the muscle and not the *moment* delivered by the muscle (e.g., Feldman 1986; Gribble et al. 1998; Günther and Ruder 2003; McIntyre and Bizzi 1993; St-Onge et al. 1997; see also *Eq. 1*). To adequately reconstruct the EP trajectory that gave rise to the movement, it is necessary to estimate the neural input to each individual muscle from the measured movement. To calculate the neural input to a muscle, it is at least required to know how active state of a muscle changes with CE length, how muscle force changes with CE contraction velocity, and what the time delays in the reflexive pathways are. Furthermore adequate reconstruction of the EP would also require knowledge of all other inputs to the system (e.g., reference velocity, co-contraction level) and feedback gains. As this information cannot be retrieved adequately in an in vivo experiment, experimental determination of the true EP trajectory is doomed to fail. This means that alternative ways of experimentally testing EP control need to be devised, for example by comparing predicted neural inputs to the muscles with recorded electromyography (EMG) or by using the claimed equifinality of EP-controlled systems (e.g., Hindler and Milner 2003).

In summary, our simulation study has shown that the stiffness estimated using a KBI model was different from the true stiffness of our musculoskeletal model. We attributed this to the neglect of higher-order components of a controlled musculoskeletal model when approximating its behavior with a second-order KBI model. Our model will clearly not give a perfect representation of some aspects of the behavior of the real musculoskeletal system, such as the force response to muscle lengthening. However, because the dynamics of the real musculoskeletal system will definitely not be less complex than that of our musculoskeletal model, it is safe to conclude that experimentally estimating stiffness using a KBI approach, as presented in the literature, will not reflect the stiffness of the true system. Apart from the problems inherent in stiffness estimation, we have shown that the KBI approach does not adequately reconstruct the EP trajectories that were used as input for the EP controller in our simulation study, and this will be all the more true for human-subject experiments.

The results of our study imply that the rejection of EP-control theories on the basis of experimentally estimating stiffness using a KBI approach and reconstructing EP trajectories, as has been done in several studies (e.g., Gomi and Kawato 1996, 1997; Wolpert and Ghahramani 2000), was unjustified. On the other hand, the finding that fast single-joint movements can be generated on the basis of EP control obviously does not imply that these movements are in fact controlled accordingly. To address whether the brain uses EP control to generate movement, EP-controlled musculoskeletal models must be extended to incorporate knowledge about the neural control systems that interact with the musculoskeletal system. In addition, it needs to be established whether EP control is feasible for mechanical systems with multiple joints and for more demanding tasks such as movements in external force fields. In our opinion, it seems too early to close the debate on whether the brain uses internal models or equilibrium point control (or both) to generate goal-directed movements.

## Acknowledgments

We thank L. Rozendaal for valuable discussions.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 by the American Physiological Society