Conclusions on motor control depend on the type of model used to represent the periphery

Within the field of motor control, there is no consensus on which kinematic and kinetic aspects of movements are planned or controlled. Perturbing goal-directed movements is a frequently used tool to answer this question. To be able to draw conclusions about motor control from kinematic responses to perturbations, a model of the periphery (i.e., the skeleton, muscle–tendon complexes, and spinal reflex circuitry) is required. The purpose of the present study was to determine to what extent such conclusions depend on the level of simplification with which the dynamical properties of the periphery are modeled. For this purpose, we simulated fast goal-directed single-joint movement with four existing types of models. We tested how three types of perturbations affected movement trajectory if motor commands remained unchanged. We found that the four types of models of the periphery showed different robustness to the perturbations, leading to different predictions on how accurate motor commands need to be, i.e., how accurate the knowledge of external conditions needs to be. This means that when interpreting kinematic responses obtained in perturbation experiments the level of error correction attributed to adaptation of motor commands depends on the type of model used to describe the periphery.


Introduction
In the field of motor control, it is controversial which kinematic and kinetic aspects of the movement are planned and/or controlled. Or what signals are sent from the brain to the periphery defined as the skeleton, muscle-tendon complexes, and spinal reflex circuitry. To investigate these signals or motor commands, a frequently used method is perturbing goal-directed arm movements with a mechanical load. (Bhushan and Shadmehr 1999;Debicki and Gribble 2004;Gottlieb 2000;Gribble and Ostry 2000;Hinder and Milner 2003;Izawa et al. 2008;Kistemaker et al. 2010;Smeets et al. 1990, etc.). These researchers used the changes in movement kinematics in response to a mechanical perturbation to conclude what aspects of the movement are controlled and what information is used to customize motor commands to the external conditions. Yet, to make these deductions on the basis of the kinematic responses observed one needs a model of the periphery.
The periphery is characterized by high-order dynamics with nonlinearities in the tendon compliance and in the intrinsic muscle properties (force-length relationship and force-velocity relationship). It is generally a good strategy to keep models as simple as possible. Should one take into account all complexity in a model for the periphery if one wants to study which variables are controlled in goal-directed arm movements? Or is a simple model adequate?
There are clear examples why an elaborate description of the intrinsic muscle properties in models of the periphery might be essential from the control of whole-body movements such as jumping. Performing a squat jump requires a pattern of joint torques that depends critically on the initial joint angles. The fact that we can perform squat jumps from various starting postures might suggest that we use precise information and internal models (Kawato 1999;Wolpert and Ghahramani 2000) to determine the required joint torques. However, it has been shown that the force-velocity relationship of the muscles compensates the effect of perturbations of starting position before a jump to a large extent (van Soest and Bobbert 1993). The robustness of their model of the periphery to this type of perturbation could lead one to conclude that muscle commands do not need to be accurately customized to starting position. A similar result has been obtained for dealing with variations in inertia in a whole-body lifting task (van der Burg et al. 2005).
Not all researchers use a model that includes an elaborate description of the intrinsic muscle properties (Sainburg et al. 1999;Scheidt et al. 2005;Shadmehr and Mussaivaldi 1994). This raises the question to what extent conclusions on motor control depend on the level of simplification with which the dynamics of the periphery are modeled. We will test four types of models of the periphery as used in previously published research on motor control to describe the elbow joint. We will investigate the responses to changes in external conditions while keeping the motor commands the same (robustness).
The simplest way to model the periphery is by neglecting all neural and muscle-tendon properties, resulting in a torque-driven model: motor commands prescribe directly the joint torques over time (e.g., Kawato 1990;Sainburg et al. 1999). The simplest way to take into account the viscoelastic properties of the muscles and tendons is to describe the periphery as a second-order time-invariant linear massspring-damper system (e.g., Scheidt et al. 2005;Shadmehr and Mussaivaldi 1994;Thoroughman and Shadmehr 1999). This second-order model is of course still a simplification of the higher order dynamics of the periphery, which might affect conclusions based on such a model (Kistemaker and Rozendaal (2011). Still other researchers model the nonlinear excitation and contraction dynamics and the tendon compliance (e.g., Nijhof and Kouwenhoven 2000;van Soest and Bobbert 1993) and the most elaborate models also incorporate the muscle spindle feedback (e.g., Kistemaker et al. 2006;Song et al. 2008).
In this study, we do not aim to evaluate which of these models best describes the periphery. This will be by definition the least simplified model if we assume that model parameters are chosen adequately. Also, we do not aim to evaluate which model leads to highest robustness against mechanical perturbation. Rather, we wish to see to what extent simplifications in the model of the periphery affect the robustness of the model's behavior, and thus to what extent robust perturbation responses may be misinterpreted to be caused by central motor control in case an inadequate model of the periphery is used. If the model's robustness to a perturbation is not altered by a simplification, then we consider this simplification to be justified. Our approach is thus opposite to the usual approach to make models more complex if the simpler model does not suffice to explain experimentally observed behavior.

Models
We used four types of models of the periphery (Fig. 1), which differed in the level of simplification with which the dynamics of the periphery are incorporated. In the most general form, we can describe the periphery in case of the elbow joint as: where T act is the net joint torque resulting from the actuator system (muscle-tendon complexes) based on the input signal (motor commands), I is the inertia of the lower arm relative to the elbow flexion/extension axis, and T ext is the external torques working on the lower arm.
As the most simplified model of the periphery, we used a skeletal model or a torque-driven model (T-model: Figs. 1a and 2a, green lines). This model did not describe the visco-elastic properties of the periphery since it ignores muscles, tendons, and muscle spindle feedback. The input for this model was joint torque over time. This means that the commanded joint torque (T com ) equal the T act in Eq. 1. We described the lower arm and hand as one segment with a segment length of 0.26 m, a segment mass of 1.65 kg, a distance of the centre of mass from the elbow joint of 0.18 m, and an inertia relative to the center of mass of 0.025 kg m 2 .
In the second most simplified model, the periphery was described as a second-order linear mass-spring-damper system (KBI-model, Figs. 1b, 2b, blue lines), as for instance used by Scheidt et al. (2005) and Shadmehr and Mussaivaldi (1994). The KBI-model strongly simplified the nonlinear visco-elastic muscle properties, the nonlinear elastic tendon properties, and the muscle spindle feedback by describing them as a linear time-invariant system with a constant stiffness coefficient (K ) and a constant damping coefficient (B). For the one-dimensional case, this resulted in where T com is the commanded joint torque over time as given by the input signal, θ is the planned joint angle over time, θ is the planned angular velocity over time, ϕ is the actual joint angle,φ is the actual angular velocity, andφ is the actual angular acceleration. Strictly speaking his model can be regarded as a forced mass-spring-damper system with the active input being planned kinematics (θ,θ) over time as well as commanded joint torque (T com ) over time.
Alternative models can be found that specify active input as either planned kinematics (De Lussanet et al. 2002) or T com  (Uno et al. 1989). As these inputs are interchangable all these models will show the same dynamic behavior to perturbations for a certain set of values for K and B. The advantage of defining active input with both the commanded joint torque and the planned kinematics is that active input can be derived by using the kinematics and the corresponding joint torques of the planned movement trajectory. Based on Scheidt and Ghez (2007) and Shadmehr and Mussaivaldi (1994), we choose K = 16 Nm/rad and B = 2.4 Nms/rad (as measured by Mussa-Ivaldi et al. 1985) and for simplicity we did not modify these coefficients for the different simulated perturbations. A similar strategy was chosen for the model parameters of the other models in this study. As our third model, we used a musculoskeletal model (STIM-model, Fig. 1c) as described by Kistemaker et al. (2006). The model incorporated the nonlinear activation dynamics, nonlinear intrinsic visco-elastic muscle properties (force-length relationship and force-velocity relationship), and the nonlinear tendon compliance (Fig. 2b, red lines). The muscle parameters were slightly modified (see Table 1) compared to Kistemaker et al. (2006) so that the maximal isometric torque-angle relationship predicted by the model fitted the curves measured by Pinter et al. (2010). The input for this model was muscle stimulation over time.
The least simplified model of the periphery used in this simulation study was an EP-controlled musculoskeletal model (EP-model: Fig. 1d) developed by Kistemaker et al. (2006). His hybrid model combines both the α-version (Bizzi and Abend 1983;Hogan 1984) and the λ-version (Feldman 1986) of EP theory. Furthermore, the model incorporated the nonlinear activation dynamics, the nonlinear intrinsic viscoelastic muscle properties, the nonlinear tendon compliance, and the time-delayed muscle spindle feedback that contributed to the visco-elastic behavior of periphery. As it has been reported that spindle feedback is inhibited at the start and end of a movement and maximal at the time to peak velocity (Shapiro et al. 2002(Shapiro et al. , 2009, the contribution of spindle feedback to the muscle stimulation was gated by means of a Gaussian function (Kistemaker et al. 2006, Eq. B2). The input signal prescribes an EP by a set of open-loop muscle stimulations that lead to equilibrium of forces at a planned joint angle and the corresponding set of muscle lengths when equilibrium in planned joint angle is reached. For each planned joint angle, multiple EP's can be calculated all leading to a different amount of stiffness in the joint angle. In line with Kistemaker et al. (2006), all EP's were calculated under the constraint that open-loop joint stiffness in the equilibrium was maximal. The planned movement was given as an input signal consisting of three EP's: the starting angle, an angle midway, and the planned angle at the endpoint. This double-step EP-trajectory was similar to the trajectory used by Kistemaker et al. (2006) and chosen this way since it has been indicated that EP does not shift instantaneously from starting to endpoint (Bizzi et al. 1982). For simplicity, feedback gains were kept at constant values of −0.4 m −1 for muscle fiber length and −0.15 s/m for muscle fiber contraction velocity in all simulations described in this study.   Table 1 Muscle parameters used in the EP-model and the STIM-model: maximum isometric force (F max ), optimum length of contractile element (l ce_opt ), slack length of series elastic element (l se_0 ), slack length of parallel elastic element (l PE_0 ), coefficient of tendon displacement method as described by Grieve et al. (1978) (a 0 , a 1e , and a 2e )

Perturbations
As shown by our definition of the periphery (Eq. 1), we can apply a mechanical perturbation by changing inertia or changing the external torque. We simulated three perturbation types: (1) a position-dependent external torque perturbation (as if moving in the gravitational field), (2) a velocity-dependent external torque perturbation (as if moving through water), and (3) an inertial perturbation (as if transporting an object from one point to the other in the horizontal plane). For all models, an input signal was obtained that would lead to a planned joint angle (θ) or a planned joint movement (θ (t)) when simulating the unperturbed situation (reference simulation). This input signal was used without modification when the perturbation was applied.
For the position-dependent external torque perturbation, we introduced gravity, as shown in Fig. 3a. We calculated steady-state joint angle given this input signal and with an additional gravitational torque (T ext ) applied at the lower arm: where m is the segment mass, g is gravitational acceleration, and d is the distance from flexion/extension axis of the elbow to the lower arm center of mass. This procedure was performed for a range of joint angles (45 • < θ < 135 • ) chosen such that gravitational torque in the perturbed condition ranged from zero (θ = 45 • ) to maximal (θ = 135 • ). We consider it suitable to use the unperturbed situation as a reference because then deviations are indicative of the contributions of the dynamics of the periphery in the responses to perturbations as used in experiments. We will refer to the difference between θ and the steady-state joint angle as the static error. For a positional task such as a goal-directed movement, this static error is in our view a crucial variable, as it indicates the perturbation-induced error in reaching the planned endpoint in the absence of additional control from the brain.
To explain the differences in static error between the four types of models, we calculated the low frequency stiffness of the STIM-model and the EP-model for 45 • ≤ θ ≤ 135 • using the methods as described by van Soest et al. (2003). In brief, this entailed linearizing the nonlinear time-invariant state space STIM/EP-models for a certain θ (linearization point). For this linearized state space model, a low frequency stiffness was estimated by imposing a constant small angular deviation from the linearization point and calculating the steady-state joint torque. For the velocity-dependent external torque perturbations and the inertial perturbations, we first simulated a 90 • elbow flexion and extension with the EP-model (unperturbed reference). For the STIM-model, T-model, and KBI-model the input signal was determined that would result in ϕ-trajectories identical to those in the unperturbed reference. For all four models, we simulated an elbow flexion and an elbow extension given this input signal in the presence of a velocity-dependent external torque (T ext ): We varied b in such a way that damping changed from −108 to +108 % compared to the damping coefficient of the KBI-model (i.e., −2.6 to 2.6 Nms/rad). A negative damping coefficient corresponds to a velocity-dependent supply of movement energy. We choose to simulate this condition as it was used to perturb arm movements in several previous studies. (Hinder and Milner 2003;Kurtzer et al. 2005;Thoroughman and Shadmehr 1999) We also simulated an inertial perturbation by using the input signal of the reference simulation and a changed moment of inertia (ranging from −50 to 125 % of the inertia in the reference simulation).
For most models, the static endpoint error is insensitive to velocity-dependent torque perturbations and inertial perturbations. To capture to what extent a perturbed movement differed from the unperturbed movement over a reasonable timeframe, we developed a measure of the model's robustness based on the simulated movements as described in phase plots. In the phase plots, we normalized position to the movement amplitude of the reference movement and angular velocity to the peak angular velocity of the reference simulation. We calculated for each time-sample in the simulation the distance in normalized state space between the reference simulation and the perturbed simulation and took the average of this distance over the first 600 ms of the simulation (Rb). Figure 3b shows the static errors induced by the gravitational torque: the difference between the actual steady-state joint angle and the angle θ . The gravitational torque increased monotonically from zero for θ = 45 • (i.e., lower arm vertical) to a maximal value for θ = 135 • according to Eq. 2. Differences in the static error among the four types of models can be explained by this relationship between gravitational torque and joint angle in combination with the differences in low frequency stiffness' of the models (Fig. 3c). For the T-model (low frequency stiffness of zero), no steady-state angle was reached if gravitational torques were exerted on the lower arm (θ > 45 • ). For the KBI-model (a linear stiffness coefficient independent of θ), the static error is approximately proportional to the size of the gravitational torques. Due to the nonlinear position dependency of the gravitational torque, the static error is not by definition proportional to the size of the gravitational torque. For the STIM-model (low frequency stiffness varies with θ) the static error increased every θ > 45 • will lead to a static error >45 • because steady-state joint angle will be reached at the upper boundary of the range of motion of the elbow joint (set at 180 • ). c The low frequency stiffness for the four types of models with gravitational torque but saturated at θ > 80 • because low frequency stiffness rapidly increased for θ > 80 • . This rapid increase in low frequency stiffness is comparable with the results of Kistemaker et al. (2007) and can be attributed to the rapid decrease in torque as found in the maximal isometric torque-angle relationship for the elbow flexors in this range of elbow angles (Pinter et al. 2010). For the EP-model, the muscle spindle feedback added to the low frequency stiffness, reducing the static error substantially. In summary, researchers who use the KBI-model or the T-model to simulate movements in the gravitational field, would come to the prediction that motor commands need to be accurately customized to the gravitational torques working on the lower arm in order to make movements that end on target. Researchers using the EP-model would come to the opposite conclusion since this model predicts that neglecting the gravitational torques would only lead to a small (<2 • ) static error.

Velocity-dependent external torque perturbation
The trajectories for velocity-dependent torque perturbation during elbow flexion simulations are shown in Fig. 4a-d for a range of perturbation sizes (−2.6 Nm/rad < b < 2.6 Nms/ rad) and for each of the four models. Results for two perturbation sizes (b = 0.6 and b = −0.6) are shown in Fig. 4e-f. These perturbation sizes are used to illustrate the difference in robustness to the velocity-dependent external torque for the four types of models ( Fig. 5; Table 2).
Again the T-model showed the least robustness to this perturbation type. No matter how small the perturbation, b < 0 always led to an oscillation with increasing amplitude and b > 0 always led to the movement ending in its starting position (see Appendix A for mathematical explanation of the latter phenomenon).
For the three models that incorporated visco-elastic muscle properties, the effect of the velocity-based perturbation was less dramatic than for the T-model. As expected, all models are less robust when velocity-dependent torque is destabilizing (b < 0) than when it is stabilizing ( Table 2). The STIM-model showed least robustness to the destabilizing perturbation and for b < −0.6 Nms/rad the intrinsic viscous muscle properties could no longer dissipate the kinetic energy injected by the velocity-dependent external torque. The result is an oscillation with increasing amplitude. The addition of muscle spindle feedback to the model of the periphe- Phase plots indicating the model's robustness to the velocitydependent perturbation for two perturbation strengths and two movement directions. Thin black curves indicate the unperturbed condition for all models (reference simulation used for normalization) ral motor system (yielding the EP-model) increased robustness to the velocity-dependent torque perturbation in both the stabilizing and destabilizing conditions ( Fig. 5; Table 2).
For b < −1.1 Nms/rad, the EP-model can no longer dissipate the kinetic energy injected by the destabilizing velocity-dependent torque. The KBI-model is also robust to the velocity-dependent torque. Only when b < −2.4 Nms/rad, the oscillations increase in amplitude. At this point, the damping coefficient of the perturbation annihilated the damping coefficient of the KBI-model.
Altogether, this means that when using destabilizing velocity-dependent external torque perturbations to investigate motor control, researchers using a T-model would predict that motor commands need to be perfectly customized to the external damping conditions in order to be able to control movements. As both sensory information as motor commands have some degree of noise, this prediction would necessitate that motor commands are adapted during the movement based on feedback, even if no perturbation is present. When using one of the three other models, the level of simplification has an important effect on interpreting data from experiments with a negative-simulated viscosity: the STIM-model predicts different responses than the other two.

Inertial perturbation
The trajectories for inertial perturbation during elbow flexion simulations are shown in Fig. 6a-d Table 2) to the inertial perturbation for the four types of models.
The movements produced by the T-model were also most affected by this type of perturbation. The three models that include damping and stiffness properties all seem to show a certain level of robustness to inertial perturbations (Fig. 6). When looking at the kinematic responses of these three models in more detail ( Fig. 7; Table 2), the model describing these properties with the least simplifications (EP-model) was least robust.
For the interpretation of inertial perturbation experiments in terms of knowledge needed to reach the target, the level of simplification in the model is not very relevant, provided that the stiffness and damping properties are incorporated in some way. If experiments would show that an inertial perturbation does not affect reaching the goal, any model that incorporates the stiffness and damping properties would lead to the conclusion that this performance is possible without knowledge of the perturbation.

Discussion and conclusions
We found that the amount of information about the external conditions that is required to control goal-directed arm movements varied with the type of model used to represent the periphery. Models that neglect the visco-elastic properties of the periphery predict that movements are not well controlled unless motor commands are perfectly customized to load conditions, whereas models that incorporate these properties predict that movements remain reasonably well controlled when motor commands are less accurately customized to load condition.

Conclusion on motor control depend on the level of simplification with which the periphery is modeled
We showed that for the model with no visco-elastic properties (T-model), all perturbations no matter how small changed kinematics substantially (no robustness). This means that this type of model predicts that motor commands need to be perfectly customized to external conditions and thus that to control movement detailed knowledge of the external conditions (inertia, gravity, and velocity-dependent external torques) is needed in addition to information on the ongoing movement to adapt motor commands during movement execution.
Researchers use this type of model to interpret kinematic responses as measured in perturbation experiments would attribute all load compensation they find to adaptation of motor commands. This affects the conclusions made on motor control. For example, if it would be found that for an unexpected 25 % change in inertia subjects show kinematic responses such as shown in Fig. 7 for the STIM-model (red curves), using a T-model to explain these data would lead to the conclusion that the adaptation of net elbow torque to the new inertia was completely due to adaptation of motor commands (for instance, by means of supraspinal feedback). In the case that a STIM-model or a KBI-model would be used, one could conclude that load compensation could occur under motor commands that are not customized to the new inertia, attributing the adaptation of net elbow torque to the dynamics of the periphery.
We found that including damping and stiffness properties in the model of the periphery improved robustness to perturbations tremendously compared to using a model that neglects them (T-model). The predictions on how accurate motor commands need to be with regard to external conditions differed also for the three types of models that included damping and stiffness properties. Therefore, when interpreting kinematic responses as measured in perturbation experiments, the amount of error correction attributed to adaptation of motor commands based on knowledge of external conditions depends on what type of model is used. This means that the conclusions one will make on motor control will depend on the type of model used to describe the dynamic behavior of the periphery. We use the position-dependent torque perturbation as an example (Fig. 2b). If in an experiment setup in which gravitational direction would be manipulated and it would be found that participants showed endpoint errors such as shown in Fig. 2 for the EP-model (grey curves), using a KBI-model to explain these data would lead to the conclusion that motor commands had been adapted to the new gravitational load, for instance, based on sensory information on gravitational direction (vestibular system, vision, kinesthetic). Using an EP-model would lead to the conclusion that all load compensation could be attributed to the dynamics of the periphery, leading to the conclusion that motor commands are not customized to gravity.
The difference in kinematic responses for the three models that include damping and stiffness properties is most prominent in the position-dependent torque (gravitational) perturbation (Fig. 3b). It is clear that the stiffness of the KBI-model has a lower value than the low frequency stiffness found for the STIM-model and the EP-model. First, is the choice of parameters responsible for the differences in static error that we found? We could increase the stiffness coefficient (K) for the KBI-model, which would increase the model's robustness to this type of perturbation. Joint stiffness measured in perturbation experiments varies between studies and higher values have been reported (de Vlugt et al. 2006;Lacquaniti et al. 1993). Yet, even with a higher stiffness, the model's response to perturbation (static error) would remain approximately proportional to the gravitational torque and therefore show a qualitatively different response than the STIM-model and the EP-model. Second, one could also argue that gravity leads (in addition to a position-dependent torque) also to a position-independent increase in torque. This might affect the behavior of the nonlinear models. We performed additional simulations in which we introduced static torques (ranging from −10 to 10 Nm), and found that this additional static torque did not affect our results in describing robustness to position-dependent torque perturbation.
The robustness to the gravitational perturbation is inversely related to the level of simplification with which the dynamical properties of the periphery are described (Fig. 3b). This is, however, not the case for the other perturbations, as is most clearly demonstrated in the velocity-dependent torque perturbation. The following order from most robust to least robust to this type of perturbation: EP-model, KBI-model, STIM-model, and T-model. Since the dynamical properties of the periphery are more simplified for the KBI-model than the STIM-model, this does not resemble the following order for the level of simplification. For the inertial perturbation, we can make a similar conclusion on following order of the robustness: KBI-model, STIM-model, EP-model, and Tmodel.
As we show that conclusions on motor control drawn from kinematic responses to perturbations depend on the type of model used to describe the dynamical properties of the periphery, one could ask which model to use in order to come to correct conclusions. This is a question that has no simple answer. If we assume that the least simplified model best describes the response of the human motor system we would suggest to use this model. Yet, we did find that for some perturbation types (Figs. 4f, 6e-f), a simplified model such as the KBI-model would describe the kinematic responses to perturbation equally well in terms of robustness compared to less simplified STIM-model or EP-model. For these perturbation types, a KBI-model would be the most adequate. The main conclusion of this study is that when using kinematic responses as found in perturbation experiment to come to conclusions on motor control, the level of detail used to describe the periphery should be considered carefully. The contributions to robustness caused by the dynamics of the periphery to the kinematic responses should be clear before drawing conclusion on central motor control.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.