Abstract
Vertebrates and arthropods are both capable of load compensation during aimed limb movements, such as reaching and grooming. We measured the kinematics and activity of individual motoneurons in loaded and unloaded leg movements in an insect. To evaluate the role of active and passive musculoskeletal properties in aiming and load compensation, we used a neuromechanical model of the femur–tibia joint that transformed measured extensor and flexor motoneuron spikes into joint kinematics. The model comprises three steps: first, an activation dynamics module that determines the time course of isometric force; second, a pair of antagonistic muscle models that determine the joint torque; and third, a forward dynamics simulation that calculates the movement of the limb. The muscles were modeled in five variants, differing in the presence or absence of force–length–velocity characteristics of the contractile element, a parallel passive elastic element, and passive joint damping. Each variant was optimized to yield the best simulation of measured behavior.
Passive muscle force and viscous joint damping were sufficient and necessary to simulate the observed movements. Elastic or damping properties of the active contractile element could not replace passive elements. Passive elastic forces were similar in magnitude to active forces caused by muscle contraction, generating substantial joint stiffness. Antagonistic muscles co-contract, although there was no motoneuronal coactivation, because of slow dynamics of muscle activation. We quantified how co-contraction simplified load compensation by demonstrating that a small variation of the motoneuronal input caused a large change in joint torque.
Introduction
During aimed limb movements, the musculoskeletal system transforms motoneuronal activity into muscle forces, which move the limb toward the goal. In an invertebrate model for goal-directed limb movements, the scratching behavior of locusts, a hindleg is aimed toward a tactile stimulus on the ipsilateral wing (Berkowitz and Laurent, 1996; Matheson, 1997; Dürr and Matheson, 2003). The targeting accuracy is maintained even if the tibia is loaded so as to increase its inertia 12-fold (Matheson and Dürr, 2003). Moreover, this increase of mass does not affect kinematic variables such as movement velocity or initial direction, indicating that the leg muscles must generate significantly more force to execute the same movement. To understand limb targeting and load compensation requires an understanding of the sensorimotor transformations that are involved in movement control. These transformations include the activation dynamics of muscles, intrinsic muscle and tendon properties, moment arms, and mechanical properties of the joints (Zajac, 1993; Pandy, 2001). Few studies integrated these factors to demonstrate their influence on control of behavior (for review, see Dickinson et al., 2000), because the complexity of the parameter space makes it difficult to identify and measure all relevant parameters in an unrestrained animal.
Neuromechanical modeling can simulate the natural movement of a limb using a defined input, and sensitivity analyses can evaluate the relative importance of measured physiological variables. Biomechanical models have been applied to examine the stabilizing mechanical properties of limbs during cockroach walking and running. In cockroach running (Full and Ahn, 1995; Jindrich and Full, 2002), a model comprising only a spring-and-damper system illustrates how the physical properties act to stabilize the animal against mechanical disturbances. The human musculoskeletal system has been similarly analyzed during locomotion (Wagner and Blickhan, 1999; Geyer et al., 2003). In frogs, the muscles stabilize wiping movements without neuronal activity (Richardson et al., 2005). Brezina et al. (2000b) transformed a neural input signal into rhythmic movements of Aplysia feeding behavior using a neuromuscular model, which was also analyzed to determine optimal neural patterns (Brezina et al., 2000a).
We used a biomechanical model to determine the role of musculoskeletal properties for the generation of unconstrained goal-directed limb movements in a prominent arthropod model organism, the locust (Burrows, 1996). We analyzed the motor behavior of hindleg scratching (see Fig. 1), during which the animal detects a target location on the wing surface and moves its hindfoot toward the stimulus position. This behavior permits a detailed analysis of the interaction of neurophysiology and biomechanics because it comprises a continuum of targeted movements of a limb with extremely well known joint and muscle anatomy (Hoyle, 1955a; Burrows and Horridge, 1974; Heitler, 1974) and completely described motoneuron innervation (Hoyle, 1955a; Sasaki and Burrows, 1998). This allowed us to model the complete neuromechanical transformation from motoneuron recordings to joint movement. Comparison of trials of the same animals with unloaded and loaded legs allowed us to evaluate the same model under different force regimens and to investigate the significance of individual model components for load compensation.
Experimental setup and example of locust scratching behavior. A, The locust was tethered with its hindleg standing on a support (black square), which determined the start posture. Stimulation at either an anterior or a posterior location on the wing (triangles) initiated a scratching movement. The hindleg performed a targeted movement in which the end of the tibia was aimed at the stimulus location, followed by a cyclical component (trajectory shown as a black line for a scratch aimed at the posterior target). In half of the trials, a load was placed at the end of the tibia (circle). A three-dimensional optical motion capture system analyzed the locations of six reflective markers (crosses) and thus measured five joint angles (B) in each video frame. Joint angles φ and ψ were 0 when the leg was in a plane parallel to the sagittal body plane. C, Synchronous electromyograms were recorded of extensor (top channel, slow extensor tibiae, SETi) and flexor (bottom channels, proximal flexor, distal flexor) motor activity. Gray shading marks extensor burst activity.
We find that a model using passive elasticity and damping is sufficient to simulate tibial movement equally well as a model with additional active properties. The long time constants of muscle activation and deactivation result in strong co-contraction of the antagonists, leading to substantial joint stiffness. We show how both of these aspects affect neural control of goal-directed limb movements and load compensation.
Materials and Methods
Experimental setup and stimulus protocol.
Experiments were performed on female locusts (Schistocerca gregaria Forskål) taken from a crowded colony at the Department of Zoology, University of Cambridge (Cambridge, UK).
Locusts were tethered above a foam ball using a fine wire noose placed around the pronotum. This allowed the animal to freely move its legs and adjust the posture (Matheson, 1997). The eyes and ocelli were blacked out using water-based black acrylic paint (Cryla, Daler-Rowney, UK). The left mesothoracic leg was amputated at the base of the femur to prevent the locust gripping the myogram wires (see below), and the stump was sealed with a wax plug (Lactona Surgident, Philadelphia, PA). A foot rest positioned posterior to the metathoracic coxal joint and ventral to the abdomen was used to maintain the same start position before each leg movement. Retroreflective circular markers of 0.5 or 1 mm diameter were cut from reflective tape (Scotchlite; 3M, St. Paul, MN) and attached to specific locations on the hindleg and thorax (Fig. 1) using clear nail lacquer. The forewing was notionally divided into five regions of equal length. The most anterior and posterior regions were used as the two stimulus sites and correspond to sites 1 and 2 (anterior) and 4 and 5 (posterior) used by Dürr and Matheson (2003). Light tactile stimulation using a fine paintbrush within one of the defined stimulus regions was used to elicit grooming behavior. Posterior and anterior stimuli were applied randomly until at least 20 scratches were recorded for each stimulus site.
In loaded condition, a weight of 2 mm diameter and 142 mg was attached to a distal tibia position, ∼2 mm proximal to the tibia–tarsus joint. The animals were filmed from a lateral view, using a color video camera (TK-C1380; JVC, London, UK) with standard PAL resolution, operated at a shutter speed of 
s. The images were taped on a VHS video recorder (NV-HS9000; Panasonic, Secaucus, NJ), displayed on a monitor, and played back for capture using a personal computer video interface card (miroVIDEO DC30 plus; Pinnacle Systems, Mountain View, CA).
We filmed N = 5 animals and used n = 80 trials (16 trials per animal, four trials for each unloaded and loaded condition as well as anterior and posterior stimulus). If more than four trials were recorded for each condition, we selected the trials with the largest number of cyclic grooming movements. Statistics were calculated using the statistics toolbox of Matlab 6.4 (MathWorks, Natick, MA). Average values are given as medians with the interquartile range (IQR).
Motion capture of scratching movements.
The recorded videos were deinterlaced, and a threshold filter was applied to obtain the marker positions in each video frame. Joint angles for each frame were computed using motion capture software (Zakotnik et al., 2004), using a three-dimensional kinematic model of the hindleg. The model contained five rotational degrees of freedom (DOF): the body–coxa joint was a ball joint with three DOF (retraction/protraction, pronation/supination, and levation/depression), the coxa-trochanter joint had one DOF (levation/depression), and so did the femur–tibia (FT) joint (extension/flexion). The tibia–tarsus joint and the tarsus were omitted in the model, because they do not contribute to aiming the tibia (Dürr and Matheson, 2003).
A disambiguation algorithm (Zakotnik and Dürr, 2005) determined the joint angle time course for each of the five joint axes. A quintic smoothing spline then used this accuracy measure to smooth the data within measurement accuracy (Zakotnik and Dürr, 2005). In any time step, the smoothed joint angles did not deviate more than 2° from the unsmoothed angles. Example joint angles for a scratching movement are shown in Figure 1.
Electrophysiology and identification of tibial motor neurons.
Muscle activity was recorded from the extensor and flexor tibiae muscles using myogram electrodes made of pairs of 20-μm-diameter polyester-insulated copper wires that were passed from the recording sites in the metathoracic leg to the pronotum, fixed at intervals to the cuticle by small drops of beeswax. The tips of the wires were stripped of their insulation to improve the recorded signal. In each animal, three recording sites were selected based on the known innervation of muscle fiber bundles and the sites of bundle insertion, which are identified by characteristic patterns of cuticular coloration.
Myogram electrodes were inserted into the distal anterior flexor tibiae muscle bundle a11, a proximal flexor tibiae bundle (p1 to a2; [nomenclature from Sasaki and Burrows (1998)]), and into the most distal extensor tibiae muscle bundle (Hoyle, 1978). Motoneurons were identified using standard criteria, including spike amplitudes, resistance reflex responses to imposed movements of the tibia, and correlation of activity with spontaneous movements (Field and Burrows, 1982; Burrows, 1995). In each animal, we recorded at least one slow and one fast motoneuron from each site, amounting to at least 6 of 11 excitatory motoneurons known to innervate the tibial muscles. Nine excitatory motoneurons innervate the flexor tibiae muscle, each one innervating only a restricted subset of muscle fibers. For example, the most proximal and distal fibers are innervated by nonoverlapping subsets of motoneurons (Sasaki and Burrows, 1998). Some fibers in the proximal bundles receive innervation from as many as seven of the nine excitatory motoneurons (Sasaki and Burrows, 1998). In addition to the nine excitatory motoneurons, the tibial flexor muscle is also innervated by two common inhibitor motoneurons (Pearson, 1973; Hale and Burrows, 1985; Wolf, 1990), which we did not record.
The metathoracic tibial extensor is 88% larger than the flexor (Bennet-Clark, 1975). Its pennate muscle fibers are innervated by only four motoneurons: the excitatory fast (FETi) and slow (SETi) extensor tibiae neurons, a common inhibitor (CI1) (Hoyle, 1955a,b; Usherwood and Grundfest, 1965), and a modulatory dorsal unpaired median neuron (Hoyle et al., 1974; Bräunig and Pflüger, 2001). Individual muscle fibers are innervated in 1 of 12 neuron combinations (Hoyle, 1978). The fibers themselves are classified as slow, intermediate, or fast based on their ultrastructure and thus contractile properties (Hoyle, 1978). Fast fibers are innervated by FETi only, whereas slow fibers are innervated by only SETi. Intermediate fibers are innervated by both FETi and SETi (Hoyle, 1978).
Electromyogram (EMG) signals were amplified (1000×) using standard extracellular amplifiers and filtered with custom-built 50 Hz notch, 500 Hz high-pass and 5 kHz low-pass filters. Signals were captured directly to computer using a Cambridge Electronic Design 1401 interface and Spike 2 software (Cambridge Electronic Design, Cambridge, UK). EMGs were displayed off-line using a software custom written in Matlab 6.5 (MathWorks). Individual spikes were detected using a combination of adjustable thresholds and manual identification based on criteria described above. Slow and intermediate flexor motoneuron spikes on a single channel were treated as a single pool. For the proximal flexor site, we classified only the spike with the largest amplitude as a fast flexor spike, although these muscle bundles should be innervated by two fast motoneurons.
Occasional crosstalk or artifacts were discounted manually by careful inspection of all recordings. SETi, FETi, and flexor motoneuron spikes were then analyzed separately for each recording site. Different slow motoneurons were recorded in both flexor channels, but they were active at approximately the same time. Only trials without fast motoneuron (FETi or fast flexor neurons) activity were used in the modeling described here, which were 61% of all recorded trials.
Muscle model and joint mechanics.
The transformation of recorded motoneuron spike trains into joint torque was modeled by including the musculoskeletal properties, the apodeme (tendon) insertion points, and the mechanics of the joint system. Each muscle was modeled as a sequence of two functions: activation dynamics and muscle contraction. The activation dynamics described the time course of isometric force given a certain spike train of motoneuron action potentials. In particular, it transformed a sequence of action potentials, α(t), into a continuous active state of the muscle, γ(t), which may be interpreted as the number of actin–myosin bridges. Function α(t) was 0, except for brief periods when an action potential occurred, when it resembled a half sine wave with a width of 1 ms, as proposed by van Zandwijk et al. (1996). Activation dynamics were modeled by two coupled overcritically damped differential equations (Eqs. 1, 2) based on Hatze (1977) and J. Zakotnik and L. Blackburn (unpublished observations). 
β(t) can be interpreted as the time course of the muscle fiber membrane potential.
A force potentiation factor c(f) was introduced to model calcium-dependent facilitation. It depended on the instant spike frequency f and increased the force for spike frequencies around 20 Hz. It was defined as a Michaelis–Menten-type equation: 
Activation dynamics therefore contained six parameters (θ1–θ4 and K1–K2), as determined by previous experiments done by Zakotnik and Blackburn (unpublished observations). These parameters were kept constant in the present study. Parameters θ1–θ4 and K1–K2 were derived from mean parameter values of a model for isometric force production on the extensor muscle when stimulated by the slow motoneuron SETi (Zakotnik, 2006). The values were θ1 = 79, θ2 = 2783, θ3 = 4919, θ4 = 78582, K1 = 1.46 × 10−2, and K2 = 3.9 × 10−4. Using these parameters, a single twitch had a rising time of 54 ms and an exponential decay time of 110 ms to half of peak twitch force. Although these parameter set was not optimized on force measurements of flexor muscle, flexor activation dynamics measured by Phillips (1980) suggest that flexor parameters should be very similar to extensor parameters.
To account for the slow muscle dynamics, muscle activation was calculated 100 ms before the scratching movement, which began when the locust first moved its tarsus off the foot rest. In addition, the active state γ in Equation 1 was scaled by 8.15 × 10−5 so that an isometric force of 20 mN corresponded to a stimulation frequency of 50 Hz, as measured by Blackburn (2005). The parameter set used for Equations 1 and 2 was also used for the flexor muscle, but the scaling factor was adjusted to account for lower measured flexor force. For both muscles, the maximum isometric force could be adjusted by the optimization procedure (see below). Activation dynamics were calculated separately for spike trains from the three recording sites, so the distal and proximal flexor recordings yielded two separate time courses, which were then averaged to obtain overall flexor active state γf.
Extensor force Fe and flexor force Ff were obtained from the muscle active states γe and γf and a passive force component arising from a parallel elastic element. The active force component was obtained by scaling the active state by the maximum isometric force F̂ and by two muscle characteristics that depend on relative muscle length l and muscle shortening velocity −v (Zajac, 1989): 
In Equation 3, function L denotes the force–length (FL) characteristic, function V is the force–velocity (FV) characteristic, and function P is the passive force. Relative extensor length le and flexor length lf were set according to the standard muscle model proposed by Zajac (1989) to be 1 at an FT joint angle of 90°, 0.5 for 0°, and 1.5 for 180°. The maximum joint angle range of the FL relationship served as an approximation for an average muscle fiber, because the true muscle fiber lengths for different joint angles are not known and may be heterogeneous in different parts of the muscle (Hoyle, 1955a). For intermediate angles, the muscle length was derived separately for extensor and flexor from the geometry of the FT joint (Heitler, 1974), as described later in Materials and Methods (Eqs. 12, 13).
The FL characteristic L reduced muscle force for extreme muscle lengths. It was a Gaussian-type function as used by Geyer et al. (2003), with adjustable width w: 
During optimization, parameter w determined the decay of active force for muscle lengths beyond the resting position. The FV characteristic V was a sigmoid function scaled to the interval [0;2] with parameter s, which controlled its slope: 
The tanh function approximated the original function described by Hill (1953) but had only one parameter controlling its slope in contrast to two lumped hyperbolas in Hill's model. The tanh function in Equation 6 behaves essentially like the pair of Hill parabolas in that a quickly contracting muscle does not produce any force, whereas a lengthening muscle could maximally double the maximum isometric force (compared with a factor of 1.8 used by Zajac, 1989). During optimization, parameter s determined the decay of force during contraction. The passive force P of the musculotendon system was implemented as a cubic parabola controlled by a scaling parameter p, depending on muscle length l (Eq. 7). The cubic parabola resembles the passive force function used by Zajac (1989): 
During optimization, parameter p was varied to set the relative contribution of passive elastic force to total muscle force production.
To calculate the joint torque τ, extensor force Fe and flexor force Ff were multiplied by their moment arms Me and Mf, which depended on the FT angle γ: 
Moment arms were given by the geometry of the joint and the attachment points of the apodemes. The relationship between joint angle and moment arm (Eqs. 9, 10) was derived from Heitler (1974), their Figures 1 and 2, and geometrical relationships described for the FT joint by Field and Coles (1994), their Appendix. The extensor moment arm Me changed approximately sinusoidally with the joint angle, whereas the flexor moment arm Mf was more complex, because the apodeme was redirected over Heitler's lump and the tibia was bent at the insertion point. The flexor moment arm was largest for flexed positions and decreased monotonically at larger joint angles: 
In this model, the flexor moment arm Mf would become negative for joint angles larger than 154°, because the apodeme remains redirected over the lump, despite the fact that this is not the case in the animal. Nevertheless, the model was appropriate for the present study because all measured FT joint angles were lower than 122°.
The geometry was also used for the calculation of relative muscle length depending on the FT angle γ: 
The viscosity of the joint system, caused by cuticular damping and passive forces of the musculoskeletal system, was modeled as a torque opposing the joint movement. Therefore, the net torque τn was given by the torque τ produced by muscles and a velocity-dependent damping with proportionality factor d: 
The net torque τn was the torque that actuates the FT joint in the forward dynamics simulation.
Dynamics model and simulation.
To derive the FT joint movement given the calculated FT torque and movement of the other joints and segments, a three-dimensional forward dynamics simulation was implemented using Newton–Euler equations (Spong and Vidyasagar, 1989). In the simulation, the FT joint was actuated by the torques obtained by the muscle model (Fig. 2) and moment arms, whereas the other joints were moved passively to reflect exactly the recorded joint angle time courses during each scratch. This allowed us to consider torques induced by passive FT joint and tibia segment movement attributable to proximal segment movement and gravity.
Scheme of the biomechanical model used to simulate locust scratching movements. Numbers indicate the order in which the components of the model were included, thus increasing its complexity in steps. The neuronal input consisted of slow motoneuron action potentials (SETi and slow flexor), which were transformed into the active state of the muscle using activation dynamics. The active state was scaled by the maximal possible isometric force and modified by passive elasticity, force–velocity, and force–length properties of the muscle. Muscle force was transformed into FT joint torque by moment arms given by the joint geometry. A joint damping component reduced the torque, depending on joint velocity. Tibial movement was calculated from joint torque using a forward dynamics algorithm. Co-contraction occurs when both antagonistic muscle models produce a nonzero torque.
The simulation was implemented using Matlab Simulink (MathWorks), which uses a variable-step Runge–Kutta integrator (ode45) to integrate joint angles from calculated joint accelerations. On a 3 GHz Pentium IV computer, the model calculated ∼0.8 simulation seconds per real-time second. Time courses of the FT joint in the simulation were sampled at 5 kHz, and the recorded movement of proximal joints was interpolated from 50 Hz to 5 kHz. In addition to the kinematic model described in the previous section, the simulation also incorporated the mass of the tibia and its inertia tensor. Segment masses and lengths were measured and set for each animal model separately. The mean value for all tibia masses was 14.9 mg (SD of 2.3 mg) and for all lengths 22.8 mm (SD of 0.63 mm). The inertia tensor was that of a homogenous cylinder with radius 0.5 mm. When modeling loaded experiments, a simulated weight of 142 mg was added to the tibia inertia tensor, using Steiner's theorem. The weight was modeled as a cylinder with a length of 2 mm and a radius of 0.5 mm placed at the distal end of the tibia and in its center.
Optimization protocol for model variants.
Up to six model parameters were optimized with respect to the least-squares error between the simulated and the recorded FT joint angle time courses of each trial. The complexity of the five model variants was increased cumulatively, as depicted in Figure 2. The simplest model variant 1 scaled maximum isometric extensor and flexor force separately; thus, it contained two parameters. This variant assumed that the isometric muscle force acted linearly on the moment arms, which transform force into joint torque. The second model variant 2 additionally comprised a viscous joint damping component (Eq. 14), in which a joint torque was subtracted from the sum of antagonistic torques (three parameters total). This variant assumed a linear transform from isometric into isotonic muscle force but included an opposing velocity-dependent frictional torque. The third model variant 3 included a passive elastic force (Eq. 7) with a common parameter for both muscles (four parameters). This variant thus took into account a passive elastic property of both muscles, rendering the musculoskeletal system into a passive spring-and-damper system. The fourth model variant 4 included an FV characteristic (Eq. 6) for both muscles with a common slope parameter (five parameters). Variant 4 therefore included an active damping property of the contractile element, which was not present in variant 3. The fifth variant 5 contained an FL characteristic (Eq. 5) that modulated the force production depending on muscle length (six parameters). Thus, variant 5 included an active-spring property of the contractile element. Other model variants were tested as controls as indicated in the corresponding sections in Results.
Optimization was performed using a bounded Levenberg–Marquardt algorithm implemented in Matlab (lsqnonlin). To prevent the algorithm from converging to local minima of the error function, the optimization was restarted with random values at least 60 times. The numerical bounds within which parameters were initialized and could vary by were determined such that the range of solutions was not limited by these bounds (Table 1). We used measured values as an estimate when available. These estimates were used to limit the range in which parameters were initialized: maximum isometric extensor force with SETi stimulation at 50 Hz was measured by Blackburn (2005) and set to F̂e = 20 mN. Because the measurements were done at a joint angle of 140°, the true value of the model parameter F̂e is likely to be larger. For the flexor muscle, the same range was assumed, although the muscle is weaker than the extensor. However, not all motoneurons were recorded; therefore, an increased isometric muscle force could compensate for a lack in recorded neural activation.
Optimization constraints for model parameters
Measurements of joint viscosity are not available for the FT joint, so we referred to an experiment on cockroach leg performed by Garcia et al. (2000), in which the cut leg was moved in a pendulum-like oscillation. We chose the initial damping parameter such that it produced a qualitatively similar damped motion of the simulated locust leg. To account for the different animals and leg mass in our study, the range of possible initial values was set to a large interval (Table 1).
Passive torques in the FT joint were measured by Burrows and Horridge (1974) for one animal, allowing us to use the joint geometry to calculate passive forces in each muscle and to fit one model cubic parabola to both the extensor and flexor muscles, giving a mean parameter value of p = 0.7.
Measurements of the FL- and FV-function parameters are not available, so their values were set to qualitatively resemble the standard muscle model proposed by Zajac (1989) and their parameters could vary across a large interval during optimization.
Results
Motor patterns driving scratching movements
Scratching behavior consisted of a targeted movement toward the stimulus position followed by one or more cyclical movements near the stimulus position (Fig. 1) (Dürr and Matheson, 2003). The targeting movement lasted 100–200 ms and consisted of a levation of the femur and extension of the tibia. The cyclical movement typically included both extensions and flexions of the tibia. During scratching, the FT joint could move from a maximum flexion of 12.9° to a maximum extension of 121.2°, with a median joint angle of 58.3°. Femur–tibia joint angles depended on the stimulus position. For the anterior stimulus position, the tibia was mostly flexed (49.2°; IQR of 26.6°), whereas posterior scratches resulted in a more extended tibia (82.4°; IQR of 32.2°). When the leg was loaded, movement patterns remained indistinguishable from unloaded movements, as judged from the direction and velocity of the initial targeting movement and the velocity, number of loops, and area covered by the cyclic part of the movement (Matheson and Dürr, 2003).
Because trials with fast motoneurons were excluded, tibial movements were driven only by SETi and slow flexor motoneuron activity. Each scratch contained three (IQR of 2) SETi bursts, with a median instantaneous spike frequency of 57 spikes per second (sp/s) (IQR of 36 sp/s), three (IQR of 4) proximal slow flexor bursts with a spike frequency of 75 sp/s (IQR of 68 sp/s), and one (IQR of 4) distal slow flexor bursts with a frequency of 66 sp/s (IQR of 69 sp/s). Although spike frequencies were higher than the rate required for maximal muscle activation [50 sp/s (Hoyle, 1978)], bursts typically contained fewer spikes than necessary for full activation of the muscle. In some cases, even single spikes occurred, which could be interleaved between bursts of an antagonistic neuron. There was no overlap between bursts of antagonistic motoneurons. In loaded trials, burst duration and spike frequency increased in extensor motoneurons, whereas both of these variables decreased in flexor motoneurons. As in unloaded trials, bursts of antagonistic motoneurons did not overlap in loaded trials (J. Zakotnik, K. L. Page, T. Matheson, and V. Dürr, unpublished observations).
Dynamics of recorded scratches
Median net FT joint torque during unloaded scratches was 1.1 μNm (IQR of 0.7 μNm), which was calculated using inverse dynamics. To compensate for gravity, almost all postures (97%) required positive torque at the FT joint, i.e., torque driving tibial extension. In loaded scratches, the mean FT torque increased 16-fold to 17.0 μNm (IQR of 12.6 μNm) to generate the movement with a mass attached to the tibia. Femur–tibia joint torque was also influenced by interaction torques caused by the motion of more basal segments, which give rise to centripetal and coriolis forces. For example, fast levation of the femur at the thoraco-coxal or coxo-trochanteral joints could cause a negative torque in the FT joint, resulting in a passive flexion of the tibia. To quantify the passive FT torque generated by body, coxa, and femur movement only, inverse dynamics were calculated in the absence of gravity and tibial acceleration. The mean absolute interaction torque for an unloaded tibia was −4.1 × 10−3 μNm (IQR of 3.2 × 10−2 μNm), which amounted to only 0.4% of the net torque. Even for loaded trials, interaction torque amounted to only 1.1% of the net torque (−0.2 μNm; IQR of 0.9 μNm).
Behavior simulation using biomechanical models
To determine the minimal model configuration sufficient for simulating scratching behavior, time courses of FT joint angle were generated using biomechanical model variants of increasing complexity (see Materials and Methods) (Fig. 2). In each new model variant, all parameters were reoptimized, regardless of whether they had been optimized in a previous model variant. One typical recorded time course is shown in Figure 3 together with time courses simulated by different versions of our biomechanical model. The error statistics for all 80 recorded and simulated trials are shown in Figure 4.
Typical recorded (dotted line) and simulated (solid line) FT joint angle time courses for different model variants, all of which included at least active isometric extensor and flexor force. The example movement was aimed at a posterior stimulus and had a loaded tibia. Top rows in each graph show the times of recorded slow motoneuron spikes (from top to bottom: SETi, slow proximal flexor, slow distal flexor). A, The simulation result with a PD component (RMSE of 29.6°). B, The time course for a model including PD and PE (RMSE of 6.1°). C, The time course including an additional force–velocity–property (RMSE of 7.0°). D, An additional force–length property (RMSE of 9.1°). In E, the PE was removed from the model variant including all other properties (RMSE of 28.4°). The small oscillations at second 1 and 2 are a simulation artifact, because the tibia moved against the physiological joint angle range constraint. In F, the PD was removed from the model with all properties (RMSE of 11.1°).
RMSE for different muscle model variants (increasing from left to right), divided into 40 unloaded (open boxes) and 40 loaded (gray boxes) trials for five animals. Boxes show the lower and upper quartile, and whiskers represent 1.5 times the interquartile range. Using only extensor and flexor forces resulted in oscillatory movements and therefore created the largest RMSE of all models (column “force”). PD reduced the error, but a PE component decreased the error further, because it returned the tibia from extreme positions to the resting position. Addition of a force–velocity (FV) and force–length (FL) property did not reduce the error.
The most simple model variant 1 using only activation-dependent extensor and flexor forces moved the tibia between extremely extended and flexed positions and, therefore, could not reproduce the recorded movements (Fig. 4). For this model, the median of root mean squared error (RMSE) for trials in both loading conditions was 30.9° (IQR of 13.3°). Addition of a passive damping (PD) component (model variant 2) decreased tibial velocity and reduced the median RMSE to 18.9° (IQR of 8.9°). However, model 2 could not generate fast movements like the cyclic component of scratching behavior (Fig. 3A). When a passive elasticity (PE) component was added (model variant 3), there was a good fit between simulated and recorded joint angle time courses (Fig. 3B). The median RMSE was reduced to 5.3° (IQR of 3.4°), which was the lowest error of all models analyzed. The error appeared to be lower when simulating an unloaded leg (4.7°) than when a load was added (6.1°), but this difference was not significant (one-way ANOVA; F(40,40) = 3.02; p = 0.09). Also, there was no significant difference between model performance for trials with anterior and posterior stimulus positions (one-way ANOVA; F(40,40) = 0.85; p = 0.34).
Incorporating FV and FL characteristics into the model (variants 4 and 5) did not further reduce the median error (FV, median of 6.3°, IQR of 4.0°; FL, median of 5.7°, IQR of 3.7°). Neither error was significantly different from the error of model variant 3 that lacked these active muscle properties (one-way ANOVAs; FV, F(80,80) = 3.27, p = 0.07; FL, F(80,80) = 2.02, p = 0.16). Example time courses for model variants 4 and 5 are shown in Figure 3, C and D. An initial extension of the tibia in Figure 3D was attributable to the scaling property of the FL characteristic, which strongly reduces flexor force during the initial phase of the trial because muscle length is very short in this flexed position. At the same time, passive elasticity of the extensor causes a positive net torque at the joint. The same effect is also seen in Figure 3F.
These results, summarized in Figure 4, showed that elastic and damping properties of the musculoskeletal system significantly improved the model performance but that additional active elasticity and damping did not improve performance. To test whether the passive properties could be replaced by the corresponding characteristic of the active contractile element, we removed either the PD or PE component from model variant 5. If PD was removed, the RMSE increased significantly (one-way ANOVA; F(80,80) = 20.57; p < 0.01) to 7.9° (IQR of 4.6°), shown in Figure 3F. It extended too rapidly from flexed positions, which resulted in a large error in the initial targeting movement and overshot at the turning points of the cyclical movement. In addition, a phase shift between recorded and simulated time course was observed in many trials (e.g., after t = 1.8 s in Fig. 3F). Without the PE (Fig. 3E), the tibia did not return from extremely flexed or extended positions in two trials and was therefore not optimized systematically for all other trials.
Parameters of optimized models
A separate set of optimized model parameters was determined for each trial, so we could evaluate the parameter variability for different loading conditions and model complexity. Because all model parameters are thought to be constants of the musculoskeletal system, they should not vary with the experimental variables stimulus position and load. Rather, the optimization procedure should find one global optimum for each animal. Furthermore, because all parameters were initialized at random before optimization, small variability in a particular optimized parameter indicates that this value is important for movement generation. Such random initialization would cause unimportant parameters to be optimized to different arbitrary values for each different movement.
The median parameter values and 10,90 percentiles for each model are shown in Table 2. These overall medians are also subdivided into separate values derived from loaded and unloaded scratches (Table 2, numbers in parentheses). On average, the parameter variation (determined by 10,90 percentiles) was 59% of the median values for all parameters and model configurations. In model variant 2 that included passive elasticity, the flexor force F̂f was optimized to the lower bound of its range in most trials (Table 2). Nevertheless, lowering the boundary did not improve the joint angle RMSE, indicating that optimization attempted to minimize flexor force in this model.
Median musculoskeletal model parameter values for each model variant
In model variant 3, which comprised both passive components (PD and PE) and had the best performance, the maximum isometric forces were 55 mN for the extensor muscle and 42 mN for the flexor muscle. Addition of the FV characteristic in model variant 4 reduced the flexor force significantly to 33 mN (one-way ANOVA; F(80,80) = 7.78; p = 0.006), suggesting that velocity-dependent force reduction affected the flexor more than the extensor. Moreover, the variability of the extensor and flexor maximum isometric force decreased significantly when the FV property was included (Kolmogorov–Smirnov test with median-centered distributions, p < 0.02). At the same time, the median joint damping coefficient d was reduced from 0.6 to 0.4 (Table 2), which indicated that the FV characteristic accounted for part of the total joint damping that was necessary to achieve a good model fit. Including both active components (FV and FL characteristics) in model variant 5 confirmed the low flexor force of the previous model variant. The apparent reduction of the FV parameter was not statistically significant.
The characteristics of the optimized model variant 5 are shown in Figure 5. In general, the variance of all model parameters indicated that particular features in the joint angle time course led to strong modification of individual parameters. The shape of the active FL and FV characteristics (Fig. 5C,D) appeared to be affected more than the passive characteristics (Fig. 5A,B). The extreme cases of the FV characteristic range from a nearly linear solution to a nearly discrete switch from factor 0 to factor 2.
Parameter distributions for the complete muscle model as used in variant 5 (solid lines, median; shadowed areas, 10,90 percentiles). A, Passive joint damping, which is proportional to joint velocity. B, Passive elasticity parabola for the musculotendon system. C, Force–velocity property of the muscle (negative velocities correspond to contracting muscle). The x-axes for the joint and muscle velocity graphs are scaled to maximal velocities observed in scratching behavior. D, Force–length property of the muscle.
In all model variants that included PD and PE (models 3–5), muscle parameters did not change significantly with load, except for maximum extensor force, which was larger for loaded trials in all model variants (Table 2) (one-way ANOVA, p < 0.01, separately for each model variant). For example, the maximum extensor force for the model variant 3 including only PD and PE was 23% larger in loaded conditions than in unloaded conditions.
Generalization capability and parameter sensitivity of the model
Having determined the parameter ranges for each model variant, it was necessary to evaluate how well each particular solution generalized to all other recorded trials. A sensitivity analysis was used to determine the extent to which variation of a particular solution affected model performance.
The ability of the model to generalize across trials within one animal was evaluated by generating a standard model with a median parameter set from all trials except for one test trial. The test trial was simulated using the standard model and compared with the recorded FT joint angle time course for the same trial. This procedure was repeated for each trial, rendering the number of evaluated trials n = 80.
The error generated by the standard model on the single test trial decreased with increasing model complexity (Fig. 6, white boxes) in a pattern that was similar to the trend for individually optimized trials (compare with Fig. 4). Median generalization errors were above 25° in model variants 1 and 2 and approximately halved with the introduction of passive damping. Accordingly, the model variant 3 with PD and PE had a median RMSE of 13.1° (IQR of 8.0°). To illustrate the performance of this model variant, the inset in Figure 6 shows a representative trial (compare with Fig. 3). The generalization error was significantly larger for loaded trials than for unloaded trials (one-way ANOVA; F(40,40) = 6.4; p = 0.01), indicating that the model did not generalize well across load conditions. This reflected the significant difference in extensor force (unloaded, 50 mN; loaded, 61 mN) (Table 2). Addition of the FV property in model variant 4 caused a nonsignificant reduction in the RMSE to 12.0° (IQR of 5.9°) and abolished the significant difference between loaded and unloaded trials (one-way ANOVA; F(40,40) = 0.3; p = 0.68). This suggests that the segmentation of damping into an active and a passive component improves generalization. An additional FL property caused a nonsignificant increase in the error to 12.2° (IQR of 6.6°). Application of a standard model to new trials of the same animal increased the overall error by 6.7° on average compared with models optimized to a single trial. As for individually optimized trials (compare with Fig. 4), there was no significant difference between generalization errors for different stimulus positions in any model variant (one-way ANOVA, all p > 0.1).
Generalization and parameter sensitivity of optimized muscle model variants. Generalization (white boxes), The RMSE for each trial was computed using a muscle model determined by mean values from all other trials from one animal. Generalization performance was best for a muscle model that included passive properties and the force–velocity property. Parameter sensitivity (gray boxes), Same error measure as in generalization test but using the extreme values of the optimized parameter ranges to illustrate sensitivity to parameter variation. Sensitivity was largest for passive damping and elasticity in model variants 2 and 3. Only in these cases were errors significantly larger than in the generalization test. Inset, The sample trajectory (RMSE of 13.0°) was generated using the average optimum of model variant 3 with passive damping and passive elasticity, which was the best model for individually optimized trials.
To assess how sensitive the models were with regard to variation of the optimized parameters, we took the extreme values produced by the optimization procedure and determined the error between each recorded trial and a simulation calculated using the extreme values. In each model variant, extreme values were set for the parameter of interest, whereas the other parameters were set to their median values. The results are shown next to the generalization errors in Figure 6 (gray boxes). Parameter variation caused significant differences only in model variants 2 (one-way ANOVA; F(80,160) = 22.4; p < 0.001) and 3 (one-way ANOVA; F(80,160) = 10.0; p = 0.0018). This highlights the importance of passive damping and passive elasticity for appropriate modeling of movement.
Role of musculoskeletal model components
To assess the relative contributions of passive and active elements of the model, the net torque computed for each trial was decomposed into the torques generated by (1) the active contractile elements of the two antagonists (isometric forces, scaled by FV and FL), (2) parallel passive elasticity of each muscles, and (3) passive damping at the joint (for model structure, see Fig. 2). This analysis was performed on model variant 3 that included PD and PE and on variant 5 that included all components. Active force production in variant 3 is equivalent to the isometric force, but variant 5 incorporates both isometric and isotonic forces. The relative contribution of passive and active components is qualitatively similar in both of these model variants, so we illustrate the results for only model variant 5 because it includes all known muscle characteristics (Zajac, 1989).
Typical torque time courses for a leg in unloaded and loaded condition are shown in Figure 7. The torques generated by the extensor and flexor muscles (Fig. 7Aii,Bii, thin lines) and by joint damping (Fig. 7Aii,Bii, dotted lines) were large compared with the effective net torque acting at the FT joint (Fig. 7Aii,Bii). This was because most of the muscle-generated torque was absorbed by the antagonistic muscle and joint damping; thus, co-contraction of antagonistic muscles causes significant stiffness of the limb. During extensions, the extensor produced 53.0% (IQR of 13.9%) of the total torque, counteracted by flexor contraction (32.6%; IQR of 21.3%) and joint damping (6.5%; IQR of 13.7%). During flexions, flexor muscle produced 49.9% (IQR of 12.6%) of the total torque, counteracted by extensor contraction (47.3%; IQR of 12.5%) and damping (4.5%; IQR of 7.1%). The mean net torque was often positive (extension) even if the tibia was flexing, because the tibia needed to be lifted against gravity.
Two sample trials (Ai, Aii, unloaded tibia; Bi, Bii, loaded tibia) with the net torque (bold solid line in Aii, Bii) decomposed into the separate torques generated by the extensor (top thin line) and flexor (bottom thin line) muscles, respectively. The bold dotted line shows the torque attributable to passive joint damping. The simulated time courses of the FT joint angle are shown in Ai and Bi. Because of co-contraction of both muscles and joint damping, the net torque in the unloaded trial is 28-fold lower than the average muscle-produced torque and thus was near the 0 axis. In the loaded trial, net torque is only 1.8-fold lower.
Contribution of active and passive musculoskeletal properties
In the absence of external forces, passive muscle force (PE) moved the tibia toward a resting position of 70.9°. Although the PE is modeled as a cubic parabola for each of the antagonistic muscles, the overall passive joint torque depended approximately linearly on the FT joint angle. In other words, the joint behaved almost like a linear spring over most of the working range (Fig. 8). This linearization is attributable to the matched action of the antagonist muscles and their moment arms, both of which are nonlinear within the working range. The spring constant of the model was −0.98 μNm/°. This means that, for extreme joint angles, passive torques may reach 54 μNm, which is sufficient to lift even a loaded tibia without any motoneuron activity.
Passive elasticity of antagonistic muscles depending on the FT joint angle. Although the passive force is modeled as a cubic parabola depending on muscle length, because of the moment arms, the resulting torque resembles a linear torque spring with respect to the joint angle. Its maximum absolute torques reach 54 μNm in the model.
The relative proportions of torque generated by active muscle contraction and passive elasticity were strongly dependent on the stimulus position and the corresponding movement (Fig. 9). For movements aimed at the anterior stimulus, flexor forces were predominantly active (Fig. 9, black bars, anterior), because the FT joint was generally working in a range that was more flexed than the resting position. In other words, to keep the joint flexed required active flexion to counter the passive tendency to return the joint to the resting position. In contrast, for movements aimed at the posterior stimulus, which required operating the joint at angles more extended than rest, it was the extensor muscle that had to generate active force (Fig. 9, black bars, posterior). For example, the passive extensor torque for anterior stimulus positions and an unloaded extending tibia was 12-fold larger than the active extensor torque (median passive, 21.8 μNm; median active, 1.8 μNm). The extensor torque was opposed by active flexor torque, which, for these anterior stimuli, was larger than the active extensor torque for both extension and flexion movements. The flexor muscle therefore actively counteracted the extensor PE in flexed positions and was active not only during flexions but also during extensions (Fig. 9). A corollary of this is that there should be little motor drive to the extensor muscle during extensions from the most flexed joint angles. This was borne out in our experiments, in which some extensions were generated without or with very little SETi activity (data not shown), much like what has been described previously by Berkowitz and Laurent (1996).
Active and passive components of torque production for all trials separated by stimulus position, loading condition (0, unloaded; L, loaded), and tibial movement direction for model variant 5. The relationship between passive and active torques strongly depends on the stimulus position. In anterior stimulus position, extensor torque is mostly generated passively (gray bars), because the tibia is likely in a flexed position. In contrast, flexor generated torques are mostly active (black bars). For posterior movements, the leg is more extended than in anterior movements; therefore, the extensor passive torque contributes less to total torque, whereas the flexor passive torque is larger than for anterior movements. Net torque (white asterisks, extensor torque minus flexor torque and damping) is mostly positive to compensate gravity.
During unloaded posterior stimulus trials and a flexing tibia, the flexor active torque was low compared with the flexor PE torque, because the tibia was more extended and therefore the flexor muscle stretched (median passive, 20.3 μNm; median active, 2.8 μNm). In this case, active extensor torque was larger than the active flexor torque in extensions and flexions alike, indicating that a large part of flexions was attributable to passive flexor force (PE) and gravity.
To determine the time course of active and passive torques during the control of extensions, we evaluated 90 extensions that had a mean velocity of >100°/s and a duration of at least 160 ms (Fig. 10). This excluded short oscillating tibial movements. During first half of the extension when the FT joint was flexed (30.0°; IQR of 16.3°), passive torque was larger than actively generated torque. After 80 ms, active torque exceeded passive torque, because the leg was in a more extended position (54.0°; IQR of 26.3°). The total torque produced during these extensions was approximately constant in the model. This shows that passive forces tend to initiate extensions.
Active extensor torque (black boxes) and passive extensor torque (gray boxes) for n = 90 extensions with a minimum mean joint velocity of 100°/s and a duration of minimally 160 ms. In the first half of an extension, passive torques were larger than actively generated torques.
Co-contraction and load compensation
To extend a loaded tibia, locusts increased net joint torque by >16-fold (1627%), as calculated by inverse dynamics. In contrast, the average extensor burst for a loaded tibia generated only ∼50% more force, which can be derived from isometric force attributable to muscle activation dynamics. To examine this discrepancy, Figure 11 shows which components of the model contribute to the torque increase during extensions.
Femur–tibia joint torque during extensions is the extensor active and passive torque (left bar of each pair) minus flexor torque and damping (right bar of each pair). In unloaded trials (gray bars), co-contraction of both muscles resulted in a median net torque of 1.1 μNm, although each of the muscles produces >25 μNm. In loaded trials (black bars), extensor torque is increased by 45.8% and flexor torque is reduced by 20.9%, which results in a 16-fold increase of net torque to 18.4 μNm.
In loaded trials, the maximum extensor force parameter was increased by 22% and the level of co-contraction decreased (Fig. 11). In unloaded trials, the extensor muscle produced 26.3 μNm (IQR of 16.1 μNm) and the flexor 25.2 μNm (IQR of 15.8 μNm). The net torque, i.e., the difference between extensor and flexor torque, that was necessary to extend the unloaded tibia was only 1.1 μNm (IQR of 0.6 μNm). Therefore, the extensor produced a torque 24-fold larger than what would have been necessary without flexor co-contraction and joint damping. The observed 1627% change in net torque to 18.4 μNm could also be achieved by a symmetrical 32.2% increase in extensor torque and 32.2% decrease in flexor torque. This illustrates how large co-contraction levels increase the dynamic range of the joint torque, because a small change in muscle force can cause a strong increase in net torque. Generally, the larger the torques produced by co-contraction of extensor and flexor (τe, τf), the less is the relative change C of antagonistic torques needed to increase net torque by a value τinc: 
As one of the muscles decreases its force, co-contraction also decreases during loading. In our experiments, SETi activity increased by 30.6% and slow flexor activity is decreased by 25.6% when the tibia was loaded.
Discussion
Motor control of aimed limb movements involves the biomechanical transformation of a motoneuronal signal into movement. By modeling this transformation for aimed scratching movements in the locust (Berkowitz and Laurent, 1996; Dürr and Matheson, 2003; Matheson and Dürr, 2003), we show that the dynamic range of joint torques is increased by co-contraction. Passive properties of the musculoskeletal system cause large joint stiffness in this invertebrate model system. Both effects facilitate load compensation.
We first evaluated the necessary model complexity (Figs. 2⇑⇑⇑–6) and then the relative importance of passive and active properties (Figs. 7⇑⇑⇑–11). The minimal model that can simulate the limb movement at a low error value includes active isometric muscle forces, viscous joint damping, and passive muscle elasticity.
Model assumptions and reliability
A very complex physiological model is likely to contain redundant parameters. Moreover, molecular details may be implemented realistically, but their role in generation of movement may be difficult to determine (Dickinson et al., 2000). In contrast, we formulate a relatively simple, mathematically tractable model based on only a few plausible assumptions. It uses only single-coefficient characteristics and the same characteristics for antagonistic muscles. Our model has allowed us to determine the significance of individual characteristics that are of widespread importance in the biomechanics literature (for example, Zajac, 1989).
The model receives its input from single motoneurons, which were recorded from the EMG of the antagonistic extensor and flexor tibiae muscles. The flexor muscle is innervated by nine motoneurons (Sasaki and Burrows, 1998), approximately half of which we recorded. Therefore, there could be additional flexor activity missing in our input signal. Because the activity of slow flexor motoneurons at both recording sites was correlated, we assume that the missing motoneurons were also active during similar time periods. All excitatory motoneurons were recorded for the extensor. Our experiments sampled much of the movement range of a locust hindleg. The comparison of unloaded and loaded trials extended the dynamic range of required muscle forces beyond the typical range of scratching behavior. In the highly specialized motor programs underlying forceful kicking, the slow motoneurons of both tibial muscles can produce higher spike rates and longer bursts than we recorded (Burrows, 1995), but we are confident to have sampled most if not all of the natural physiological working range for scratching.
Each leg muscle also receives input from 1 or 2 inhibitors (CI), which were not recorded. The CI neurons receive sensory input from tarsal receptors (Usherwood and Runion, 1970), which do not have any contact during scratching. The CI neurons are active during searching movements in the mesothoracic leg (Wolf, 1990) and decrease relaxation time of the muscles. In contrast, Usherwood and Runion (1970) have shown that the CI neurons may not be active in the metathoracic extensor muscle during slow walking. If the CI neurons were active in scratching, the time course of muscle activation should have decayed faster than in the activation model used. The level of co-contraction would have been reduced but probably not abolished. During load compensation, reduced co-contraction increased net joint torque (Fig. 11). Moreover, optimization of the model to loaded trials results in a 22% increase of extensor force with no change in flexor force (Table 2). However, because the maximum muscle forces are constant in the animal, net torque could be increased by a faster relaxation of the flexor muscle. Thus, the model predicts an increase of CI activity in loaded trials.
In all model variants, the recorded neural signal was transformed into muscle force using a Hill-based muscle model with a parallel passive elasticity and a force–length–velocity property (Zajac, 1989). All model variants lacked a serial elastic element because we assumed the apodeme to be stiff. This is supported by measurements in the cockroach showing that apodemes are 40 times stiffer than a vertebrate tendon (Ahn and Full, 2002).
Within the force–length characteristic, muscle fiber length was assumed to change by at most 50% of its resting length at extreme joint positions (Zajac, 1989). If the range for the muscle length deviated from this assumption, this would decrease the PE force. Therefore, the optimized PE parameter could be biased by the muscle length range. For example, this bias could explain why the optimized PE values are 2.5-fold larger than the value estimates made by Burrows and Horridge (1974).
Because pennation angle of muscle fibers is negligible for stick insect extensor muscle (Guschlbauer et al., 2005), a muscle similar to that in our study, it was omitted in the model. Muscle activation dynamics were assumed to be the same for extensor and flexor. The time course of activation was measured for the slow extensor fibers. Because only trials with slow motoneuron activity were used, we reason that flexor activation and deactivation times were in the same range as those of the extensor.
The overall low error in the predictions of our model suggests that additional biomechanical properties that could influence force generation do not play a significant role. Even if the parameter estimates used in the model deviated from real physiological values, the model would elucidate the relative contribution of the components in movement generation.
Role of passive forces in limb movement
In arthropods, passive forces may be strong enough to carry the weight of the animal (Yox et al., 1982). In the absence of neural activity, the PE of the model acts to move the FT joint to a resting position of ∼70°. This is similar to the angle measured by Burrows and Horridge (1974). To deflect the joint from this resting position, an active torque must be generated. Therefore, the control of the movement is strongly dependent on the joint angle. For example, if the tibia is flexed beyond resting position for anterior stimuli, the extensor PE generates a torque sufficient to resume the resting position in a purely passive manner. Thus, to generate a sequence of extensions and flexions in the cyclic grooming response to anterior stimuli, only flexor activation needs to be modulated. The nervous system could therefore use the level of muscle activation to control the joint angle rather than the joint torque, although muscle activation is primarily associated with muscle force. In the simplest case, an actively produced joint torque could depend linearly on the joint angle. Figure 12 suggests that this is the case for a wide range of joint angles. Locusts were more likely to use SETi at femoro-tibial angles greater than the passive resting angle and more likely to use slow flexors at angles less than rest. It appears that nonlinearities of passive elasticity and moment arms combine to yield a linear torque, permitting the animal to rely less on motor activation at favorable angles.
Measured probability for an SETi spike (black bars) or a slow flexor tibiae motoneuron spike (gray bars) defined by the ratio between the number of spikes and the number of recorded FT joint angles in each bin. The probability for an SETi spike increased for increasing joint angles, whereas the probability for a slow flexor spike decreased for increasing joint angles.
According to our model calculations, passive forces are in the same range as actively generated forces and can thus explain tibial extensions observed without extensor activity (Berkowitz and Laurent, 1996) and our data. A similar model for passive forces could also be applied for wiping behavior of a denervated frog, for which targeting accuracy does not deteriorate when a mechanical disturbance is applied (Kargo and Rome, 2002; Richardson et al., 2005).
A second passive mechanism contributing to joint stiffness is the viscous damping included into our model. Garcia et al. (2000) have shown that internal viscosity in the joint structure and surrounding tissues considerably damps oscillations of a cockroach tibia, but they did not examine damping during natural movements. The damping property in our model could account for a linear relationship between spike frequency and joint velocity as observed in a cockroach leg (Watson and Ritzmann, 1998).
Co-contraction of antagonists
We did not observe neural co-activation of antagonistic muscles, i.e., antagonistic motoneuron bursts did not overlap, even if the additional load caused prolonged burst duration and increased spike frequency in extensor motoneurons (Fig. 1) (Zakotnik, Page, Matheson, and Dürr, unpublished observations). Nevertheless, there is a substantial amount of co-contraction of muscles because of the slow activation and deactivation of muscle fibers (Fig. 7). Slow time constants of the stick insect extensor and flexor muscle fibers have also been reported by Bässler and Stein (1996).
Co-contraction is energetically expensive, because muscle forces are absorbed by the antagonist and do not contribute to movement generation. In humans, it increases joint stiffness (Hogan, 1984; Carter et al., 1993), reduces the effect of interaction torques from the movements of other joints (Gribble and Ostry, 1999), and increases accuracy of targeted arm movements (Hogan, 1980; Gribble et al., 2003). In cockroaches, joint stiffness attributable to co-contraction is a key parameter for running speed and adaption to different surface compliances (Full and Farley, 2000) and stabilizes the impact when the leg touches ground at the beginning of the stance movement (Ahn and Full, 2002). Matheson and Dürr (2003) hypothesized that co-contraction could play a role in load compensation during locust scratching. Our model shows that co-contraction of antagonists enables the muscles to increase net joint torque with only small changes in neural activation (Fig. 11). Thus, a 43% change in extensor muscle torque is sufficient to increase the joint torque 16-fold and to move the loaded tibia against gravity in all postures. In contrast, without co-contraction, the extensor would need to generate 16-fold more torque.
Our model shows that biomechanical properties of the musculoskeletal system strongly affect the motor control of aimed limb movements. The large passive joint stiffness in the model predicts that the joint angle is a relevant parameter for the neural control of aimed leg movements. Antagonist muscle co-contraction attributable to long muscle time constants facilitates substantial load compensation with little change of neural activation.
Footnotes
-
This work was supported by Biotechnology and Biological Sciences Research Council Advanced Research Fellowship AF13046 (TM) and a studentship of the Studienstiftung des deutschen Volkes (J.Z.). We thank Keri L. Page for help with data acquisition and valuable discussions.
-
T. Matheson's present address: Department of Biology, University of Leicester, University Road, Leicester LE1 7RH, UK.
- Correspondence should be addressed to Dr. Volker Dürr, Department of Biological Cybernetics, Faculty of Biology, University of Bielefeld, P.O. Box 10 01 31, 33501 Bielefeld, Germany. Email: volker.duerr{at}uni-bielefeld.de
