Abstract

Spinal circuits may organize trajectories using pattern generators and synergies. In frogs, prior work supports fixed-duration pulses of fixed composition synergies, forming primitives. In wiping behaviors, spinal frogs adjust their motor activity according to the starting limb position and generate fairly straight and accurate isochronous trajectories across the workspace. To test whether a compact description using primitives modulated by proprioceptive feedback could reproduce such trajectory formation, we built a biomechanical model based on physiological data. We recorded from hindlimb muscle spindles to evaluate possible proprioceptive input. As movement was initiated, early skeletofusimotor activity enhanced many muscle spindles firing rates. Before movement began, a rapid estimate of the limb position from simple combinations of spindle rates was possible. Three primitives were used in the model with muscle compositions based on those observed in frogs. Our simulations showed that simple gain and phase shifts of primitives based on published feedback mechanisms could generate accurate isochronous trajectories and motor patterns that matched those observed. Although on-line feedback effects were omitted from the model after movement onset, our primitive-based model reproduced the wiping behavior across a range of starting positions. Without modifications from proprioceptive feedback, the model behaviors missed the target in a manner similar to that in deafferented frogs. These data show how early proprioception might be used to make a simple estimate initial limb state and to implicitly plan a movement using observed spinal motor primitives. Simulations showed that choice of synergy composition played a role in this simplicity. To generate froglike trajectories, a hip flexor synergy without sartorius required motor patterns with more proprioceptive knee flexor control than did patterns built with a more natural synergy including sartorius. Such synergy choices and control strategies may simplify the circuitry required for reflex trajectory construction and adaptation.

For animals without prolonged parental care there may be strong evolutionary selection pressure on rapidly elaborating useful movement. In lower vertebrates, sensing and control options may also be more constrained than those in mammals. The proprioception of frogs is determined entirely by limb state and motor outflow—there is no separated gamma intrafusal control, there is no corticospinal projection. Conceivably, the spinal cord and other motor areas may need some evolutionarily predetermined organization and capabilities—a collection of primitives or bootstrap for rapid movement development (see Giszter et al. 2001, 2007; Ting and McKay 2007).

A basic issue in producing targeted limb movements (including wiping) is that the limb has position- and velocity-dependent properties. Muscle forces depend on muscle length and velocity of shortening (Gordon et al. 1966), moment arms depend on joint angle (Pandy 2000), and the inertial load depends on limb configuration (Mussa-Ivaldi et al. 1985). Thus limb mechanics play a central role and must be modeled as accurately as is feasible. A fixed motor pattern initiated from different starting positions is likely to fail to reach a specific target from certain of those positions, due to limb nonlinear properties, although it is known that equifinal fixed patterns might exist or be discovered after deafferentation and training (Bizzi et al. 1984). Indeed, we have shown that in spinalized frogs removal of proprioceptive feedback led to target misses in reflex wiping and resulted in curved movement paths from certain positions (Kargo and Giszter 2000b).

Manipulation of proprioceptive inputs can cause significant amplitude and phase effects in primitives, but these effects preserve the pulse and muscle composition structure (Kargo and Giszter 2008). However, it is not known whether these modifications are necessary or competent to reproduce the observed limb trajectories. During wiping in spinal frogs, manipulation of proprioceptive inputs with muscle vibration causes significant amplitude and phase effects in three different individual primitives: vibration of curarized iliofibularis muscle provided independent amplitude scaling of hip and knee flexor pulses and time shifting of hip extensor pulses (Kargo and Giszter 2008). We used a detailed biomechanical model, based on Kargo and Rome (2002) and Kargo et al. (2002), to test the roles of these mechanisms. We hypothesized that if the frog spinal cord uses these several simple proprioceptive feedback mechanisms to set up trajectories, it would incorporate information about the starting limb position during unperturbed wiping at the earliest possible time and would likely use simpler computational strategies or heuristics (see Giszter and Kargo 2000; Giszter et al. 1989; Ting and McKay 2007). For example, earlier work in frogs suggests that limited simplifying heuristics may support wiping kinematic transforms for potentially complex targeting on the back. It was shown that the spinal cord of the frog is not effective if challenged to move beyond these simplifications (Giszter et al. 1989). Here we test the competency of the individual adjustment mechanisms that were identified in Kargo and Giszter (2008) and their combination to organize trajectories. We also recorded some of the proprioceptive feedback available to the frog and examined its competence to set up trajectories from different starting positions early in the wipe using these mechanisms. We sought the simplest and earliest transformations from sensory input to motor output that were competent to replicate kinematic, kinetic, and motor pattern behavior across positions. Taken together, our results show that, using appropriately chosen primitives and the simple reflex adjustments identified previously, the model limb can produce successful accurate targeting of the hindlimb from a large range of starting positions. This was possible based only on the initial bursts of muscle spindle feedback from two muscles. In contrast the model limb without such adjustments replicates deafferented frog behaviors.

METHODS

Our goal was to test the competence of fixed duration and fixed composition synergy bursts (or force primitives) to execute wipe trajectories when adjusted by one, or several, physiologically identified mechanisms, using a detailed biomechanical model. We used afferent and motor activity recording to parameterize the model.

Vertebral laminae 6–10 were removed to expose the lumbar cord and dorsal roots 7–9 on the ipsilateral side. The frog was placed on a molded stand that supported the body and hindlimbs in the horizontal plane. The pelvis and vertebral column were immobilized with custom-made clamps. Dorsal roots 8 and 9 were placed on miniature beeswax-coated platforms. Minutien (insect) pins were used to pin the fascial sheaths onto the platform (see Hoyle and Burrows 1973). This helped reduce movement-related loss of units. Small strips of tissue paper were dipped in 5% protease solution and placed on the roots for 5 min to soften the fascial sheath. Tungsten electrodes (10° taper, 5–8 MΩ; A-M Systems, Everett, WA) were used to penetrate the roots and record units. Small bipolar hook electrodes were lowered into the canal for stimulation of the ipsilateral ventral roots. Muscle spindle units were identified by several criteria (as reported in detail in Giszter and Kargo 2002). Spindles had increased firing rates with palpation of the muscle belly and imposed limb movements that lengthened the expected muscle of origin. Stimulation of ventral roots 7, 8, or 9 (ipsilateral side) produced a distinct pause in firing, when the resting firing rate was high (2–10 Hz), e.g., when the muscle was stretched. Ventral roots were stimulated between 0.1 and 0.3 V using 1-ms biphasic pulses delivered at 50 Hz. If pauses in firing rate occurred, we concluded that either the skeletofusimotor fibers (i.e., beta-motoneurons) innervating the recorded spindle were in a different, nonstimulated root or skeletofusimotor activation was insufficient to counter the effect of muscle shortening. Last, small levels of intramuscular stimulation (0.1 to 0.3 V, 2-ms biphasic pulses, 40 Hz) were delivered through the implanted electrodes to activate the presumptive muscle of origin. Spindles showed both pause and burst behavior to intramuscular stimulation and no response to stimulation of adjacent muscles.

To record muscle spindle firing patterns during wiping, we evoked the wiping response in 17 frogs by electrically stimulating the contralateral hindlimb skin at the dorsum of the ankle via bipolar leads (2- to 3-mm separation, 600-ms train, 3–8 V, 2-ms pulses, 30–40 Hz). Electromyographic (EMG) signals were led to a differential amplifier (A-M Systems) for wide band-pass filtering (100 Hz and 10 kHz) and amplified with a gain of 1K. The recording and ground leads for single-unit recordings were led to a differential amplifier where the signals were band-pass filtered (cutoffs: 100 Hz and 10 kHz) and amplified with a gain of 10K. EMGs and single-unit recordings were A-D converted using a Das16 converter and Mega FIFO memory storage (Computer Boards), sampled at 10 kHz, and stored and further filtered off-line on computer using custom-written software.

Hindlimb movements during wiping were recorded by a video camera, with a shutter speed of 1/1,000 s at a rate of 60 fields/s. The camera was placed above the frog, perpendicular to the horizontal plane of the wiping movement. A light-emitting diode (LED), triggered by the data collection program, was used to synchronize video with EMG and single-unit recordings to within <17 ms. In eight frogs, we attached an infrared LED to the ankle to more accurately measure the EMG onset to movement onset delay. An Optotrak system (Northern Digital, Waterloo, Canada) recorded the position of the LED (resolution of 100 μm) at 1,000 Hz. We found the mean EMG to movement onset delay to be 50 ms (±8 ms SD) in the eight frogs.

EMG, kinematic, and single-unit data were analyzed in Matlab (The MathWorks, Natick, MA). EMGs were rectified and filtered as described in previous studies (Giszter and Kargo 2000; Kargo and Giszter 2000a,b). The approximate joint center of the hip, knee, and ankle were digitized from video. The horizontal plane curvature of the ankle trajectory and proximity to the stimulus at the target location were determined, using the maximum deviation from the straight line path to target and a normalized path length λcurv in the two-dimensional (2D) projection of the path. This was the ratio of actual path length to the straight line Euclidean distance, as used in previous studies (Kargo and Giszter 2000b). In simulation we tested the normalized path length λcurv of the horizontal plane projection of the path and the minimum Euclidean distance λdist from target achieved by 350 ms following lift-off from the substrate

λcurv=∑i=2n(xi−xi−1)2+(yi−yi−1)2(xn−x1)2+(yn−y1)2(1)

where x and y are Cartesian coordinates, along the path that was divided into n segments. To compare with experimental data, the horizontal planar projection λcurv for simulations was used.

Instantaneous firing rates of muscle spindles were computed and plotted at the time of each spike by dividing by the previous interspike interval.

The dynamic equations of motion for the model were derived with SD/Fast (Symbolic Dynamics, Mountain View, CA). Forward dynamic simulations were produced by Dynamics Pipeline (Musculographics, Santa Rosa, CA). Simulated movements were driven by the forces generated by the 16 Hill-type actuators. For the simulations here these actuators were combined into three muscle groups or synergies, the memberships of which were identified previously by both perturbation and signal analysis methods (Giszter and Kargo 2000; Hart and Giszter 2004; Kargo and Giszter 2000b). Balance of muscles within the group or synergy were determined as described in the following text. Each muscle within a synergy received a weighted version of the same excitation signal (i.e., common drive), in a fixed pulse to form a force primitive.

The aim of this study was to examine how a large number of muscles (16) are coordinated in a simple way as three primitives to produce targeted hindlimb movements rather than to test how secondary or tertiary properties of any one muscle affects performance. Therefore simpler muscle models were implemented as a first step (see the following text). Clearly, more sophisticated muscle models would be necessary for more refined performance analysis (e.g., see Cheng et al. 2000; Sandercock and Heckman 1997). However, in our case, we used a computationally efficient, three-state Hill model. Thus we did not include a muscle series elastic element (therefore no crossbridge stiffness), no history-dependent terms, no activation-dependent shifting of the force-length and force-velocity relationships, and no motor units with different mechanical properties. The implications of not including such mechanical complexity are raised in the discussion.

Since Rana pipiens and Rana catesbaiana are considered geometrically and anatomically similar (Emerson 1978; Lombard and Abbot 1906) and since wiping in these frogs is also similar (Fookson et al. 1980; Giszter et al. 1989; Schotland and Rymer 1993; Sergio and Ostry 1993), we were confident that the EMG, muscle spindle, and kinematic data that we measured in Rana catesbaiana could be used to drive, control, and calibrate the wiping simulations, which were performed using the Rana pipiens model, for the purposes presented here.

The simulations were implemented in OpenSIMM or in Matlab-Simulink using software developed by Dr. Rahman Davoodi (MMS Software, Bioengineering Department, University of Southern California). Simulations were run to generate data in 1-ms time steps, with variable integration steps in the 1-ms interval. In the Matlab environment, the pelvis was fixed to the ground in the same position relative to gravity as experimental frogs. The limb began in the simulation positioned in the horizontal plane. Gravity was in effect throughout the simulation. Rather than simulating multiple ground contacts in the resting leg, which would be unwieldy and complex, we balanced torques due to gravity for the first 50 ms with opposing torques. These torques were applied about the hindlimb joints for the initial 50 ms of the simulation, to stabilize the chosen posture to prevent any movement arising either from force transients in the hindlimb actuators at the onset of the simulation or from gravitational forces. These forces approximately simulated the support of the substrate and static friction on this substrate at the outset of motion. At 50 ms the support surface was effectively removed and the limb had to be self-supporting and controlled in the gravitational field. This matched the measured mean delay between EMG onset and movement onset in real frogs (see earlier text). Simulations were subsequently run for 350 ms, sufficient time for a spinal frog to complete a brisk initial trajectory in the hindlimb “placing” phase (Kargo and Giszter 2000b).

The muscle forces, lengths, and activations and the joint angles, velocities, and accelerations were saved in Matlab for each simulation run. Both isometric and free-limb runs of the model were performed. The model information was then used to calculate the endpoint force at the ankle in isometric trials (i.e., when the limb joints were restrained) and the ankle trajectory during free-limb trials. Endpoint force was calculated by determining the Jacobian matrix at the restrained limb configuration and by multiplying the inverted transpose of the Jacobian matrix by a vector of static joint torques (see Kargo et al. 2002; Tsai 1999). Joint torques were determined by multiplying muscle moment arms by muscle forces. Ankle trajectories in free-limb trials were calculated using a standard forward kinematics algorithm (Denavit–Hartenberg notation) to transform the joint angles into an ankle position (see Kargo et al. 2002).

The curvature of the ankle trajectory and proximity to the target location (assumed to be midway between the hips and 8 mm caudal to the pelvis) were determined as described earlier and in previous studies (Kargo and Giszter 2000b).

Our model parameterization at this location to construct the model followed three steps. First, we adjusted the neural activations of muscles in primitives to capture our experimental data for primitives' force orientations under isometric conditions across the workspace. To find a unique set of ratios of muscle activation for each synergy, we used a procedure that minimized activations while matching force direction at a single central location in the limb workspace. We began with an exploration of the muscles observed in the synergies in the frog, but allowed muscles to drop out if the force direction was matched and the overall simulated field approximately matched the observed fields. We explored HF: ILe, ILi, ILf, GL, SA, TA KF: STd, STv, Ilf, SA; HE: GR, SM, AM, AV, STd, STv, SA, PL. We tested 10,000 different muscle balances in each group, selecting those that matched the target force directions of the synergies. We then selected among these the minimum activation sets that satisfied the force direction criteria. This procedure avoided an exploration of the full 5 DOF workspace and 16 DOF muscle space, and restricted search to a manageable region. However, it could (and did) lead to different synergies and some muscles could be dropped (e.g., SA in the HF synergy).

Second, we adjusted amplitudes and phasing of the three primitives in wiping to capture the time course of isometric force recorded experimentally in wiping at the single chosen posture.

Third, for simulation of free limb trajectories, we uniformly reduced the amplitude of the entire pattern from the second step to replicate recorded kinematics from the isometric hold posture. This uniform scaling of about 0.9 was not guaranteed to succeed, but in practice it generated close fits to experimental trajectories. However, when we allowed nonuniform scaling, the best possible scaling to fit experimental results departed slightly from uniformity. Best results were obtained if the HE and KF primitives were scaled by 0.90 (AHE = 0.90AHE isometric and AKF = 0.90AKF isometric) and the HF primitive by 0.85 (AHF = 0.85AHF isometric). This procedure provided a simulation of free-limb wiping from the single chosen posture. Subsequently, adjustments would be required based on available proprioception at other configurations unless simulating a deafferented or fixed response animal.

The individual muscle synergy pulses of fixed durations used for the first step produced viscoelastic force fields of the following form

F(r,r˙,t)=∑iNAi⋅a(t+τi)⋅Φi(r,r˙)(2)

where each A is an amplitude scaling factor, a(t) is a time-course function and Φi(r, ṙ) represents the invariant structure of the force-field primitive, which is a function of both position r and velocity. In general, earlier testing of this framework was done under isometric conditions where ṙ = 0 and can be omitted, where in that case F(r, t) represents the isometric forces that are produced at the ankle over time t and at limb position r. The isometric endpoint forces were composed of the sum of force-field primitives ϕ(r) through time. These forces will be referred to as “endpoint” forces to distinguish them from muscle contractile forces. The endpoint forces associated with the targeted portion of wiping appear to be composed of the sum of three primitives ϕ1–3(r) (Kargo and Giszter 2008). An isometric force-field primitive is defined to have the property that at a given limb position, the forces have a fixed direction over time. Force direction and peak magnitude A depend on the isometric limb position, but the time course of forces a(t) does not vary significantly across positions. During wiping, the three force primitives are activated at different latencies τ1–τ3.

Figure 1A shows the basic muscle model structure and Fig. 1B the control structure used to simulate wiping. c(t) was distributed to the muscles forming each primitive. c(t) was then scaled by Ci for each muscle. Figure 5B shows the isometric force patterns produced by an experimental frog (dashed lines) and by the model frog (solid lines) when using the above parameters and modeling structure. Both the direction and magnitude of forces during the simulated wiping response closely matched those produced by the real frog after adjustments of primitives.

Simulating hindlimb wiping with a detailed model of the frog pelvis–hindlimb complex. A: the 13 hindlimb muscles forming the model are shown as red lines. Colored arrows mark the direction of the 3 force primitives that comprise the isometric wiping response: KF (knee flexor primitive), light purple; HE (hip extensor primitive), green; HF (hip flexor primitive), dark purple. B: the framework used to simulate wiping (left to right): each primitive had a time-course generator, which output a normalized waveform (peak = 1.0) at time τ. The variable A scaled this waveform, which was then distributed to each of the muscles comprising the primitive. Each muscle had a muscle-specific variable C that scaled the excitation waveform. The muscles generate contractile forces MF that are transmitted through the limb to produce an isometric endpoint force (at one position) or force field FF (when forces are measured across a range of positions). Normalized force fields produced by each primitive are shown in the far right. When the model limb is freed to move, MFs drive the motion of the model. MF values are in turn regulated by limb motion (e.g., force–velocity and force–length properties of muscle and stress–strain properties of in series connective tissue). In our model, sensory feedback from muscles potentially regulate τ, A, or C.

Simulation details

Dynamic systems equations for wiping simulations

q¨=M−1(q)[RM⋅FM(q,q˙,lM,aM)+G(q)+V(q,q˙)](3)l˙M=f(lM,q,q˙)(4)

where q is a vector of the generalized coordinates (three hip angles and flexion-extension angles at the knee, ankle, and tarsometatarsal joints); lM and aM are vectors of muscle length and activation, respectively; M(q) and RM are the system mass and moment arm matrices; FM(q, q̇, lM, aM) is a vector of muscle forces; G(q) and V(q, q̇) are vectors of gravity and motion-dependent terms.

Muscle forces.

The force that each muscle generated FM(q, q̇, lM, aM) was determined by scaling generic muscle-tendon properties with five muscle-specific parameters: PO, peak tetanic force; lOM, optimal fiber length; α, pennation angle at optimal fiber length; lOT, tendon slack length; and εOT, tendon strain when tendon force PT = PO (see Zajac 1989 for details). Maximum shortening velocity was the same for each muscle (vO = 10.35 s−1). Lutz and Rome (1996) measured this value for semimembranosus fascicles, which were composed mainly of type 1 fibers, the fastest twitch fibers. Most frog hindlimb muscles examined have a majority of type 1 fibers that comprise 60–95% of the muscle mass (Lutz et al. 1998). Inactivated muscle was assumed to generate a passive force PPE at lengths lM > lOM and activated muscle was assumed to generate an additional, contractile force PCE. We assumed the isometric force–length relation for submaximally activated states was a scaled version of that for maximal activation—i.e., PisoCE[lM, aM(t)] = PfaCE(lM)·aM(t), where PfaCE is the force–length curve described for maximally activated [aM(t) = 1] frog muscle (Gordon et al. 1966). We assumed the contractile force to exhibit a force–velocity relationship PCE(vCE) in all muscles that was similar to that described for frog sartorius (Edman 1979). Contractile force could therefore be described by the following

PCE=[PfaCE(lM)·aM(t)/PO]·PCE(vCE)(5)

Tendon (in series connective tissue) was assumed to be elastic and have the same exponential stress (PT)–strain(εT) relation as that of the frog plantarus tendon (Trestik and Lieber 1992). Pennation angle α was assumed to be constant, resulting in the following: PT = PM·cos α, lMT = lT + lM·cos α, and vMT = vT + vM/cos α. The dynamic equation for contraction dynamics is described by the following (see Kargo and Rome 2002; Zajac 1989)

dPTdt=f[PT,lMT,vMT,aM(t)](6)

and

a˙M={[u−aM][c1u+(c1⋅⋅⋅c2)T]u≥aM[u−aM]c2u<aM}(7)

where c1 and c2 are 1/muscle activation time constant (13 ms) and 1/deactivation time constant (50 ms) (Szentesi et al. 1997) and u is the neural drive to the muscle.

Simulation of a force-field primitive and associated synergy.

Individual muscles produce torques leading to endpoint forces. These individual muscle force contributions are related to the endpoint forces associated with a synergy and its primitive by the following equations in isometric conditions (velocity terms are dropped)

ϕ(r,t)=∑iNBia(t)θi(r)(8)

where the measured forces ϕ(r, t) are the sum of the endpoint forces θi(r) produced by N different muscles. Each muscle produces endpoint forces θi(r), where θi(r) is a vector comprised of the x, y, and z components of the endpoint force normalized to the peak force magnitude for that muscle over all test configurations at peak activation. Muscles produce endpoint forces with peak magnitude Bi and normalized time course a(t). If a(t) differed among the coactivated muscles, the endpoint forces would not have a constant direction and would appear to rotate over time (Giszter et al. 1993). Thus the defining feature of a force-field primitive is that a pulsed activation of a collection of coactive and covarying muscles (a fixed synergy and burst) underlies the basic force response.

In the isometric limb, individual muscles produce contractile forces MF that are related to the endpoint forces by the following

Θ(r,t)=(JT)−1·T(r,t)(9)T(r,t)=MF[Lf,c(t),C]·r¯(10)

where the endpoint forces recorded at the ankle Θ(r, t) depend on the static hip and knee torques T(r, t) produced by muscle contraction and on the inverted transpose of the Jacobian matrix (JT)−1. JT describes the effects of limb configuration and segment lengths on the local differential properties of the limb linkage (Tsai 1999). The magnitudes of the joint torques depend on the isometric contractile force MF and the muscle's moment arms at the hip and knee (r̄). Finally, MF in the isometric limb is a function mainly of the muscle fiber length Lf at position r (relative to an optimal fiber length for force production) and the peak neural activation C of the muscle. c(t) is the normalized (neural driven) activation time course of the muscle, which varies in amplitude between 0 and 1 in the time interval.

From Eq. 8, Bi represents the peak endpoint force that each muscle produces. The set of Bi determines the direction and magnitude of the force primitive. Bi depends on the associated Ci (peak neural activation or weighting) of each muscle in the synergy (see the appendix in Giszter et al. 1993Tresch et al. 1999). However, the transformation from ϕ to a set of Bi and associated Ci for a given field structure does not have a unique solution. Further, it is difficult to uniquely predict a set of Ci for simulation purposes from EMG measurements. Thus to determine an approximate set of Ci for each synergy we did the following: 10,000 sets of randomly generated Ci (from 0 = no activation to 1 = maximal activation) were used to scale the tonic activation of those muscles previously identified physiologically as forming each synergy. The direction and magnitude of the endpoint forces were then calculated. The sets of Ci that reproduced the correct 3D force directions [i.e., those that satisfied a cos (dir [Fx Fy Fz]Experimental·dir [Fx Fy Fz]Model) <5°, where·represents the dot product, and dir is the normalized direction vector] and generated the correct force magnitude (±0.05 N) at the test location were determined. As expected, for each primitive, there were multiple synergies (sets of Ci) that satisfied these criteria. From these, we then used a static optimization, selecting the synergies that both minimized the total activation levels (i.e., ∑ Ci) and that produced forces across a set of tested limb positions that approximated the force fields measured experimentally (e.g., see the force fields in rightmost column of Fig. 1B). This produced the set of synergies used in subsequent steps. Later in the simulation we reevaluated the HF simulation synergy and relaxed the minimum activation criterion in the initial selection to better match details of observed spinal synergies.

The second step in simulating wiping was to match the isometric force pattern produced during the wiping response. This was performed at the same test position in which the sets of Ci were determined. The following Gaussian function was shown to adequately represent the time course of isometric forces produced by each primitive (see Kargo and Giszter 2000a)

a(t)=exp[−(t−ρ)2/(2σ2)](11)

where ρ is the time of peak force, t is time in ms, and σ is 40 (σ determines the width of the Gaussian, roughly 2.35-fold σ at half-peak). The rectified and filtered EMG bursts associated with each force primitive can be approximated by a similar function, i.e., c(t) is nearly the same as a(t), although time shifted. The time of peak activation (ρ) for each group of muscles were selected from published data for typical frogs [i.e., ∼70% of all frogs studied (Kargo and Giszter 2000a,b): ρKF = 90 ms, ρHF = 125 ms, ρHE = 195 ms]. Onset times were ρ − 90 (τKF = 0 ms, τHF = 35 ms, τHE = 105 ms). The other 30% of frogs exhibit a nearly synchronous burst of activity in both limb flexors and extensors (Kargo and Giszter 2000b). In this study, we focus on the typical wiping motor pattern in which primitives are activated sequentially and in a partly overlapping manner.

The third step in simulating wiping was to match the free-limb kinematics of experimental frogs from the same test position as used in the isometric force patterns. The simulated constraint at the ankle of the model was released and the same motor pattern that was used to generate the isometric response was activated. The joint angles and ankle trajectory were determined. In practice, although uniform scaling of the pattern generated good fits, we found that the best activations to match the trajectory of experimental frogs departed slightly from uniform scaling of the pattern. Best results were obtained if the HE and KF primitives was scaled by 0.90 (AHE = 0.90AHE isometric and AKF = 0.90AKF isometric) and the HF primitive by 0.85 (AHF = 0.85AHF isometric). This reduction in EMG amplitude in free-limb motion is generally consistent with the observed differences in EMG amplitude between isometric and free-limb conditions in frogs (Giszter et al. 1993; Kargo and Giszter 2000b) [Fig. 6D shows an example of the hip and knee angles and ankle velocity for the real frog (dashed lines) and the model (solid lines)]. The trajectory of the model closely matched that of the real frog up to the time of target limb contact after which the trajectories deviated. In this study, we were interested solely in the targeted portion of wiping and the model captured the isometric and dynamic characteristics of this portion from a single test position.

Tuning of model parameters.

After initial tuning to capture a single trajectory, subsequent tuning used the trajectory features λdist and λcurv as optimization criteria to adjust motor pattern features. This was done in a highly constrained fashion based on physiological data. Kargo and Giszter (2008) reported that iliofibularis spindle vibration caused specific effects in spinal frogs, independently modulating the strength of hip flexor and knee flexor activation busts, and the timing of hip extensor activation bursts. Each burst was modulated independently but all had similar durations, consistent with the decomposition of the pattern into sets of independent pulses as described in Hart and Giszter (2004) and Kargo and Giszter (2000). Semitendinosus vibration could also modulate these bursts in an analogous but slightly different fashion (unpublished data; Kargo 2000). The work of Cheung et al. (2005) and Richardson et al. (2005) suggest limited on-line correction of the ongoing trajectory without cutaneous triggering. Spindle recordings (presented in results) suggest that early in premovement muscle bursts the beta-spindle system can estimate limb state (also see Weber et al. 2007, re mammals). We thus limited our analysis to initial setup of the trajectory using the identified mechanisms and early feedback patterns. Our analysis tests the competence of these data-based descriptions to account for trajectory setup and subsequent execution. We do not intend to suggest other mechanisms are not operating during execution.

Model feedback adjustment parameters were optimized using the nonlinear least-squares optimization algorithm in the Matlab optimization toolbox (the Levenberg–Marquardt method called by m-file lsqnonlin.m), using trajectory derived variables λcurv and λdist.

RESULTS

We tested a simple control framework based on the observations of wiping trajectory controls and modularity in Hart and Giszter (2004) and Kargo and Giszter (2000, 2001, 2008). We examined whether using three pulsed primitives and the simple proprioceptive regulations described in their data could adequately capture the organization of force and kinematic patterns observed during the wiping reflex, from different positions across the workspace. We tested whether experimentally observed gain and phase shifts in the primitives (Kargo and Giszter 2008) could separately and, then in combination, appropriately modulate muscle activity to perform the observed limb trajectories based on proprioceptive feedback. Accordingly, we first examined firing patterns in the beta spindle system of the frog.

Proprioceptive information about the initial hindlimb position during wiping

Although the direct effects of muscle vibration of selected muscles and the effects of deafferentation are described (Kargo and Giszter 2000b, 2008), the pattern of spindle firing in wipes in frogs has not been published. Beta-spindle firing patterns are likely to constrain the control possibilities of the frog. To examine the pure beta-spindle patterns of the frog we recorded from hindlimb muscle spindles during nearly 1,000 wiping trials in 17 different frogs. Spindle firing patterns were examined in three time periods: resting, premovement (EMG onset up to movement onset), and during wiping limb motion. We found that the resting firing rates of all recorded muscle spindles were relatively low and ranged from 0 to 18 Hz across the configuration space of the hindlimb (65 muscle spindles from 12 hindlimb muscles; 10 ILf, 4 STv, 5 STd, 5 GR, 4 SM, 8 CR, 3 TFL, 5 GL, 6 PL, 4 ILi, 2 SA, 9 TA). In most spindles (48/65), resting firing rates averaged over a 30-s period were negligible (<1 Hz) at certain regions in the hindlimb's workspace. Resting rates were never >20 Hz. Low firing rates were likely due to reduced postural tone in the spinal frog and the absence of an independent gamma-fusimotor system. The very low resting discharges were probably insufficient to provide adequate information about the starting position unless integrated over many seconds. Because of this, we examined whether firing rates within the premovement period as beta-drive–engaged spindles carried primarily information about the starting limb position.

Only select muscle spindle populations fired strongly during the premovement period. Figure 2, B–E shows the firing patterns of four different muscle spindles during the time periods of interest (dark gray area, prewipe; light gray, premovement; white, movement). The instantaneous firing rates are shown (bottom row) in addition to the EMG activity of the muscle of origin. Figure 2, B and C shows that ILf and ST muscle spindles both fired bursts of action potentials during the premovement period. Because no obvious movement occurred yet, these bursts were due mainly to beta activation of intrafusal endings. The peak firing rates for all recorded ST and ILf spindles approached 80–130 Hz during the premovement period, depending on where the limb was located as the frog initiated wiping. Spindles from muscles activated in later phases of wiping (Fig. 2D, GR spindle; Fig. 2E, CR spindle) did not fire during the premovement period, but rather as the muscles of origin contributed to limb motion or were passively lengthened.

Muscle spindles from semitendinosus (ST) and iliofibularis (ILf) fire bursts of action potentials during the premovement phase of wiping. A: the typical kinematic pattern observed during hindlimb wiping. Only the wiping limb kinematics are shown (thigh and calf segments), from motion onset (O) to target limb contact (TC) and 7 frames after target contact (i.e., limb extension phase). Frames shown at 16.67-ms intervals. B: the rectified electromyograms (EMGs) of ILf and gracilis (GR), the discharge pattern of an ILf muscle spindle, and the instantaneous firing rate of the spindle are shown from top to bottom. The dark gray area is the resting period, light gray area is the premovement period (from EMG onset to movement onset), and white area is movement period. The dashed vertical line marks the time of target limb contact. C: the rectified EMGs of ST and GR, the discharge pattern of an ST muscle spindle, and the instantaneous firing rate of the muscle spindle are shown from top to bottom. D and E: muscle spindles from muscles not activated in the initial EMG phase did not exhibit bursts of action potentials in the premovement period; GR and cruralis (CR) spindles, respectively.

Compared with resting firing rates, the premovement firing rates of ST and ILf spindles carried more information about the initial limb position (Fig. 3). Contour plots in the top row of Fig. 4A show how mean resting firing rates of these spindles (computed over a 10-s period) varied when the ankle was placed at different workspace positions. The contour plots in the bottom row show how mean firing rates varied during the premovement period at the same limb positions. For each contour plot, the x–y values represent the position of the ankle relative to the hip (black dot represents the hip; gray oval represents the target location). Firing rates were recorded at 15 different positions and smoothly interpolated across the workspace of the limb. The premovement firing rates at one position represent the mean of four to five trials. Dark red represents a 60-Hz change from the minimum firing rate. Minimum firing rates were 0 Hz for resting ST and ILf muscle spindles and 68 and 78 Hz during the premovement for ILf and ST spindles, respectively. Dark blue represents positions at which minimum firing rates were recorded. The premovement firing rates of ST spindles were modulated over a 60-Hz range across the workspace and never dropped to negligible rates. Compare this to only an 18-Hz range in the resting limb and near-zero firing rates across certain areas of the workspace. The premovement firing rates of ILf spindles also varied considerably, over a 50-Hz range, compared with only a 7-Hz range in the resting limb.

The mean frequency of the initial burst in iliofibularis (ILf) muscle spindles varied with initial limb position. A: rasters of the spike times of ILf muscle spindle for 7 wiping trials from the same starting position (top) and the averaged firing rate are shown (bottom). B: the mean frequency of ILf muscle spindle firing was determined over 100-ms period starting at the onset of the biceps femoris (BI) EMG (top). The area under the ILf EMG envelope was determined over this 100-ms period (bottom). The mean frequency of the initial burst is shown at 5 different starting positions (hip angle = 0° for all; knee angles = −130 [most flexed]; −110, −90, −70, and −50° [most extended]).

The premovement firing rates of ST and ILf muscle spindles provide information about the initial limb position. A, top row: resting firing rates of ST (1st column) and ILf spindles (2nd column) at different starting positions across the limb's workspace (total of 15 positions; see text). Colored contours represent the change in mean frequency (over a 10-s period) relative to a minimum frequency. The minimum frequency (dark blue) for both muscle spindles was about 0 Hz. Black dots represent the hip location, the gray dot marks the target location, and bold lines represent the thigh and calf at one position. The bottom row shows premovement firing rates of ST and ILf spindles across the workspace (4–5 wiping trials per position). Minimum firing rates (dark blue) were about 60 Hz for the ILf muscle spindle and 70 Hz for the ST muscle spindle. There were larger variations in the firing rates of both spindles during the premovement period. B: the difference between the premovement firing rates of ST and ILf spindles provided a coarse representation of the orientation of the limb axis (axis defined as a straight line connecting the hip and ankle). C: the sum of the premovement firing rates ST and ILf spindles provided a coarse representation of the length of the limb axis.

The time window required to extract the same resolution of position information from ST and ILf spindles was thus clearly less during the premovement compared with that of the resting period. If need be, a rapid measure of instantaneous firing rate could be computed within a 6- to 15-ms interval during the premovement period as opposed to a 50- to 1,000 ms (interspike) interval during the resting period.

There was an interaction of muscle length and beta drive contributing to the rates in this window. We recorded from ILf and ST muscle spindles at one of five limb positions in which the knee flexor muscles were successively lengthened by incrementally increasing knee angle (130, 110, 90, 70, and 50°). We then evoked 7–10 wiping trials at each position. We found that (mean) spike frequency of the premovement burst was linearly related to muscle length at these positions (see Fig. 3). The magnitude of the (rectified and integrated) homonymous EMG at these locations was also linearly related to muscle length.

We found that the information provided by the premovement firing of ST and ILf spindles could be used to estimate the initial hindlimb position in a polar coordinate frame; specifically, the orientation and length of the hindlimb axis within the primary plane of motion (limb–axis refers to a line connecting the hip and ankle). We found the difference between ST and ILf firing rates provided a coarse description of limb orientation (see Fig. 4B) and the sum of ST and ILf firing rates provided a coarse description of limb–axis length (see Fig. 4C). Thus the essential information to extract a coarse representation of the initial limb position is provided by only two spindle populations (ST and ILf) and can be extracted within a 40- to 60-ms window at the onset of wiping. Combining feedback from other spindle populations will likely lead to significantly finer representations of limb position (e.g., see Weber et al. 2007). However, in the present study we used these data to examine whether ST and ILf information (and a coarse representation) were sufficient to set up appropriate, position-dependent motor patterns.

Fixed-pattern simulations at a single starting position

We first tested force patterns and then kinematic patterns beginning at a single starting position (see methods, Simulating hindlimb wiping with a primitive-based framework) to see whether the pattern could be generalized. We simulated hindlimb wiping using the three steps described in methods. First, we found an appropriate muscle synergy for each of the three force primitives.

The three primitives have been termed the knee flexor (KF), hip flexor (HF), and hip extensor (HE) primitives (Kargo and Giszter 2000a,b). When the hindlimb of experimental frogs is in the posture shown in Fig. 2A (hip angle θH = 0°, knee angle θK = −80°), the KF, HF, and HE primitives generate mean force directions of 190, 105, and 255°, respectively, in the horizontal plane (represented by the colored arrows; Figs. 1A and 6, A and B). From a rear view of the frog, the KF, HF, and HE primitives generate mean force directions of 28, 19, and 25°, respectively, in the frontal plane. These mean directions were determined from 10 frogs in which primitives were selectively modulated in amplitude or time of onset or were selectively deleted from the behavior (see Giszter and Kargo 2000; Kargo and Giszter 2000a,b).

Three muscle synergies account for roughly 90% of the EMG variance during the targeted portion of wiping (Hart and Giszter 2004), wipe deletion (Giszter et al. 2000), and corrections (Kargo and Giszter 2000a) and these synergies are closely associated with the three force primitives (Kargo and Giszter 2008). The muscles STv, STd, Ilf, and SA show covarying activity during the initial portion of wiping and are most often associated with the KF force primitive. The muscles ILi, ILe, GL, SA, TA, and ILf show covarying activity and are closely associated with the HF primitive (Kargo and Giszter 2000a). Finally, the muscles GR, AM, AV, SM, PL, SA, STv, and STd show covarying activity and are closely associated with the HE primitive (Giszter and Kargo 2000). For our model, we needed to find the appropriate ratios of muscle activation within these three synergies to reproduce the force primitives.

The synergies we found and evaluated in motor patterns were the minimum activation synergies that roughly captured the peak force-field patterns of experimentally observed synergies statically across the workspace (see methods). In the present study, the set of Ci for KF was 0.55·STd, 0.42·STv, 0.58·Ilf, 0.13·SA; for HE was 0.42·GR, 0.43·SM, 0.26·AM, 0.18·AV, 0.30·STd, 0.28·STv, 0.43·SA, 0.14·PL; and for HF was 0.5·ILe, 0.42·ILi, 0.18·ILf, 0.15·GL, 0.14·TA. The HF pattern discovered by this means dropped SA activation completely, although this absence was never observed experimentally. We proceeded to use this reduced set first, expecting a failure of the synergy. However, this did not occur (see following text). We then subsequently restored SA to the synergy in a final analysis.

The amplitudes and timing of the three primitives were next determined, to minimize squared error (i.e., deviation) from the trajectory of observed experimental force development measured in isometric conditions. The solid lines in Fig. 5B represent the direction and magnitude of the forces generated by a typical frog during a wiping response at a single limb position (hip angle = 0°; knee angle = 80°). The dashed lines in Fig. 5B represent simulated forces in which three force primitives were summed through time. Different aspects of these primitives are shown in Figs. 5C and 6, A and B: ϕ1(r) = [−0.872 −0.154 0.464], τ1 = 0, A1 = 0.62 N; ϕ2(r) = [−0.258 0.962 0.090], τ2 = 50 ms, A2 = 0.35 N; ϕ3(r) = [−0.257 −0.958 0.120], τ3 = 150 ms, A3 = 1.7 N; a1(t) = a2(t) = a3(t) = a(t), as described in Eq. 11. The finding that wiping forces can be adequately described as the sum of three force primitives creates a starting point from which to model this behavior in a more realistic manner (see following text).

Forces generated by the wiping limb during the targeted portion of the hindlimb to hindlimb wiping reflex can be described as the sum (through time) of 3 force primitives. A: experimental setup. The vertebral column and pelvis were fixed, the target limb was stimulated to evoke wiping in the contralateral limb, and isometric forces were recorded at the ankle of the wiping limb. The wiping limb was positioned at different start positions (see right for joint angle convention). B: isometric force direction (top) and magnitude (bottom) are shown during wiping in an experimental frog (solid lines). Modeled forces reproduce the experimental force pattern (dashed lines). C: modeled forces were generated by summing 3 force primitives (P1–P3), which were activated at different times (top); each primitive had a Gaussian time course, an onset time, a peak magnitude, and a constant force direction (bottom).

Simulating wiping forces and kinematics with the primitive framework and model frog. A: model structure. B: activation of muscles was constrained to conform to a set of synchronous synergies/primitives of fixed endpoint forces. C: the isometric force pattern produced by the model frog (solid lines) closely matched the force pattern recorded experimentally (dotted lines). D: after making minor adjustments to the isometric motor pattern (see text), the model frog also reproduced the free limb kinematics of the experimental frog. The top row shows hip and knee angles. The bottom row shows ankle velocity. Dashed line marks the time of target limb contact in the real frog. The gray area (PM) represents the 40-ms premovement period between EMG onset and motion onset that is observed in real frogs.

Our simulations of force development at a single position could thus potentially be parameterized in an open-loop fashion. Primitive amplitude and timings were fitted at this position to reproduce the force pattern (Fig. 6C) and then also fitted to the free-limb kinematics (Fig. 6D) of experimental frog wipes at this position without altering timing. The strength of activation of synergies and forces at this position needed to be reduced to capture the free-limb trajectory with the highest fidelity. A uniform scaling of about 0.9 gave a reasonable fit. However, the best possible scaling to fit experimental results departed slightly from uniform scaling of the three pulses. Best results were obtained if the HE and KF primitives were scaled by 0.90 (AHE = 0.90AHE isometric and AKF = 0.90AKF isometric) and the HF primitive by 0.85 (AHF = 0.85AHF isometric). The difference of isometric and free-limb scalings suggested a role for feedback in scaling responses.

We next sought to test a range of starting positions for free-limb motion to examine how well the pattern generalized. We intended to test whether a single pattern (after the scaling described) was sufficient to drive trajectories across the workspace. If this simplest adjustment failed we would examine how sensory feedback across positions might be used to tune the amplitude and onset times of the primitives forming the wiping motor pattern in the model.

Simulating wiping: tuning of primitives based on initial limb configuration

We tested whether an unadjusted motor pattern, which was successful at reaching the target from a single test position (see methods, Simulating hindlimb wiping with a primitive-based framework), could also reach the target from other starting positions. We first examined a range of positions at different starting orientations of the hindlimb. The length of the limb-axis from hip to ankle was kept constant (θK = −75° at each position) and the orientation was varied by changing only the hip angle (position 1, θH = −70°; position 2, θH = −45°; position 3, θH = −20°; position 4, θH = 5°; position 5, θH = 30°). Figure 7A shows ankle trajectories from these five starting orientations (300 ms of each trajectory is shown). The two black circles represent the locations of the right and left hip joints and the gray circle represents the target location. Figure 7A shows that the limb missed its target when the same motor pattern, which was successful from the test position (position 4), was used at other limb orientations.

AHF and τHE had to be tuned at different starting positions for the model to reach the target location. A: forward dynamic simulations were performed using the same motor pattern at 5 different starting orientations. At these positions, the length of the limb axis (from hip to ankle) was constant, but the orientation of the limb axis was sequentially changed. Under these conditions, the ankle trajectory (black lines) did not reach the same target location (gray dot; black dots mark the hip joints). B: AHF (amplitude of HF primitive) was varied from 0 to 2.0 in 0.2 increments at each position. The ankle still did not reach the target from positions 1 and 2. C: τHE (onset time of HE primitive) was varied from 75 to 135 ms in 10-ms increments at each position. The ankle did not reach the target from positions 1 and 2. D: AHF (amplitude of HF primitive) and τHE (onset time of HE primitive) were optimally tuned to minimize path curvature and distance to target. The best run trials are shown in which the ankle came closest to the target and had the smallest path curvature. E: wiping simulations were run from different limb–axis lengths in which limb orientation was constant. Position 3 was the same as in A–D. The limb did not reach the target when the same motor pattern, which was successful at position 3, was used at position 8. F: simple scaling of AKF resulted in small improvements from position 8 (dotted line). However, summation of a feedback component with AKF (0.30) was fully successful at reaching the target (solid lines; see text). G: ranges of variation for target proximity λdist and path curvature λcurv for the HF amplitude scaling alone, for positions 1–5. H: ranges of variation for target proximity λdist and path curvature λcurv for the HE timing variation alone, for positions 1–5.

A single pattern failed to organize trajectories of low curvature and trajectory error across the workspace. This matches the results of Kargo and Giszter (2000b): deafferented frogs show similar kinematic deficits. Accordingly, seeking a minimal adjustment, we next tested whether one of the mechanisms identified in Kargo and Giszter (2008) could adjust the test motor pattern so the target could be reached from each starting orientation or all were needed. In keeping with the physiological results with vibration only amplitude or timing, not durations of bursts, could be adjusted.

The first strategy used was to modify the peak amplitude of the HF primitive (AHF). This adjustment is observed in spinalized frogs experimentally. Figure 7B shows ankle trajectories when AHF was varied from 0 to 2.0 and in 0.2 increments at each position (10 wipes/position). When AHF was 2.0, ILe and ILi activation was near maximal. The limb did not reach the target from all starting positions when AHF was varied alone.

We tested a second strategy in which the onset of the HE primitive (τHE) was modified from 75 to 135 ms in 10-ms increments (7 wipes/position). This is the normal latency range observed across different starting positions in real frogs and an adjustment that can be controlled with muscle vibration (Kargo and Giszter 2008). The limb also did not reach the target from all positions using the strategy in which τHE was varied alone (Fig. 7C).

Our third strategy used both adjustments of strategies 1 and 2 together, as seen experimentally. In this strategy both AHF and τHE were varied simultaneously to minimize both distance to the target and ankle path curvature. We used a nonlinear least-squares optimization algorithm in the Matlab optimization toolbox (the Levenberg–Marquardt method called by m-file lsqnonlin.m) to find the appropriate input vector [AHF and τHE] that minimized the following multiobjective function

λ=wdist·λdist2+wcurv·λcurv2(12)

where λdist is the minimum distance between the ankle and the target zone during the wipe, λcurv is the ankle path curvature, wdist is a weighting factor (0.5), and wcurv = 1.0 to more equally weight the two variables. The convergence to a minimum was assessed by the algorithm. Because the optimization algorithm might converge to local minima (Neptune 1999), we started all simulations from three to five different starting points in the state space of [AHF and τHE]. We then determined whether the algorithm converged on a similar solution set from the different starting positions. Figure 7D shows ankle trajectories for “optimal” runs from each position. The parameter values that produced the best runs were the following: position 1, AHF = 1.72, τHE = 135 ms; position 2, AHF = 1.51, τHE = 122 ms; position 3, AHF = 1.00, τHE = 115 ms; position 4, AHF = 0.52, τHE = 90 ms; position 5, AHF = 0.28, τHE = 81 ms. At rostral workspace positions where the hip was flexed (position 5), AHF and τHE were smallest. At caudal workspace positions (position 1) where the hip was extended, AHF and τHE were largest.

We measured the minimum distance to the target and path curvature with each of the three strategies tested and compared these values to those measured experimentally. Table 1 shows the ranges of minimum distance between the ankle path and the target (min and max values) obtained when using the same motor pattern from each position (No), AHF modulation alone (HF), τHE modulation alone (HE), and for the best-run trials (B; simultaneous variation of AHF and τHE). The range of path curvature values (λcurv) obtained for the same conditions is also shown. The model strategy that produced best runs resulted in ankle trajectories that had a range of accuracy and curvature values similar to those produced by experimental frogs with intact feedback. However, feedback was used in a brief window only early in the wipe. The other single adjustments and the fixed patterns produced less accurate and more curved trajectories comparable to those produced by deafferented frogs.

Simulated wiping: tuning of primitives based on initial length of limb-axis

We next examined whether keeping the motor pattern fixed as the starting position was moved along one orientation axis was sufficient to reach the target. The limb started at different positions along a fixed orientation axis. Figure 7E shows the added starting positions (6–8) along the axis connecting the hip and the foot at position 3. The joint angles at these positions were as follows: position 3, θK = −75° and θH = 5°; position 6, θK = −95° and θH = 12°; position 7, θK = −105° and θH = 33°; position 8, θK = 33° and θH = −30°.

Figure 7E shows the ankle trajectories when the same motor pattern, which was successful at position 3, was used to drive limb motion from position 6 to position 8 (bold lines). This motor pattern was again not sufficient. This failure was at the most extended position, where the knee–flexor moment arms for ST, ILf, and SA were one half the moment arm values at more flexed limb positions (see Kargo and Rome 2002). An obvious solution to this problem is to enhance knee flexor muscle activation at extended knee angles, which is in fact observed experimentally (Kargo and Giszter 2000b).

We again tested the minimal strategies needed to enhance knee flexor activation and improve target success for increased polar coordinate distances. First, we scaled AKF at the extended position, thus keeping the ratios of ST, ILf, and SA activation constant. However, even with an value of 1.8 for AKF, causing maximal activation of ILf (1.0) and ST (0.99) at position 8, ankle trajectories still did not reach the target (thin, solid line in Fig. 7F, dotted lines). In the second strategy we added a uniform amount of activation to each KF muscle [FBKF·a(t), balance 1:1:1:1] in addition to the prior balance of activation levels [AKF·a(t) balance, 0.55:0.42:0.58:0.13] for the muscles STv, STd, ILf, and SA, with the latter fixed. This strategy of adding a new uniform balance of muscles FBKF modulated based on position to a fixed value of AKF was successful at reaching the target from the most extreme position (see solid lines in Fig. 7F). We found the modulation was exponentially related to polar distance (see Relationship between muscle spindle data and tuning parameters across starting positions).

We measured minimum distance to the target and path curvature at the different starting positions and with each strategy (no modulation, No; scaling of the knee flexor primitive alone, KF; summation of an optimal “feedback component” with the knee flexor primitive, B [best, meaning the value for FBKF that minimized both distance to target and path curvature, i.e., the best value]). We compared these values to experimental values. Table 1 shows the ranges of target misses and path curvature values for each strategy. The strategy of varying only FBKF was clearly most effective at the tested axis lengths and produced trajectories that had a range of target distance and path curvature values similar to those produced by experimental frogs with feedback intact.

In the previous sections, we showed that a combined adjustment strategy (combined modulation of AHF and τHE) dealt well with variations in the starting limb orientation and a second strategy (adding uniform excitation to muscles forming the KF primitive) dealt well with variations in the starting limb length. The two sets of tested positions covered two small areas of the limb's workspace. We next combined these two types of controls. Our algorithm applied in the model transforms 1) ST and ILf muscle lengths to spindle firing rates and 2) spindle firing rates to adjustments of the patterns of primitive activation. We then tested this algorithm across a large range of starting positions.

The first step of the algorithm was to transform starting muscle length to spindle firing rate. Data shown in Figs. 3 and 4 support the idea that premovement firing rates were linearly related to starting muscle length. Thus without modeling the detailed physiology or mechanics of this transformation, we simply represented firing rates as the linear functions shown in Fig. 8, A and B (data show the mean premovement firing rates ±1SD for an ILf muscle spindle and an ST muscle spindle respectively across the tested workspace [see Fig. 5]; four to five trials evoked at each limb position).

Simple use of proprioception can effectively regulate the wiping motor pattern based on starting position. A: firing rates of ILf muscle spindles during the premovement period were related in a linear manner to ILf muscle length. Firing rate data taken from Fig. 5 and muscle lengths were determined by positioning the model in the same 15 starting positions at which firing rates were recorded. B: firing rates of ST muscle spindles during the premovement period were related in a linear manner to ST muscle length. C: the difference between premovement firing rates of ST and ILf muscle spindles was related in a linear manner to the AHF and τHE values that produced the best trajectories at positions 1–5 (i.e., where limb axis orientation was varied: black dots, τHE; white dots, AHF). D: the summation of premovement firing rates of ST and ILf spindles was exponentially related to a feedback component FBKF, which was distributed more uniformly to the KF muscles (ST, ILf, and SA) than the initial fixed pattern. This produced the best trajectories at positions 3, 6, 7, and 8 (i.e., where limb axis length was varied). Equations in A–D were implemented in the wiping control and simulations were run from a large range of starting positions.

The second step of the algorithm was to combine ILf and ST feedback and transform this feedback to appropriate patterns of primitive activation. This transformation consisted of two parallel steps. First, we found that the difference between the premovement firing rates of ST and ILf spindles recorded experimentally was linearly related to limb orientations and to the values of AHF and τHE that in the model that produced accurate and straight reaches from the different tested limb orientations, e.g., positions 1–5 (see Fig. 8C). Second, we found that the sum of the premovement firing rates of ST and ILf spindles was linearly related to the tested limb–axis lengths, and the sum could be used simply to determine the values for FBKF that in the model produced accurate and straight reaches (e.g., positions 3 and 6–8; see Fig. 8D). For the FBKF strategy, the proprioception recruited SA as strongly as the other muscles and this new balance avoided saturation of ILf and ST activation. The −0.2 negative value for this strategy shown in Fig. 8D at low spindle rates indicates strong modulation of feedback across all the limb positions, including the central position at which the fixed pattern was initially developed. This strategy effectively required altered the ratios of activation of the muscles in the knee flexor group recruited due to feedback and the original synergy as an offset. This necessary addition mirrors the physiological results of muscle vibration seen in Kargo and Giszter (2008), in which SA and ST activation were regulated independently by ILf vibration (see discussion).

We tested whether the combined adjustment structure using the spindle information as shown in Fig. 8 could generate accurate and straight wiping trajectories at starting positions spanning the hindlimb's entire workspace.

The combined model in Fig. 8 can be summarized as three adjustments, of the amplitude of the hip flexor primitive, the hip extensor primitive timing, and amplitude of the knee flexor synergy, together with a fixed knee flexor primitive parameterized as follows

where ST and Ilf are muscle spindle firing rates in Hz, AHF is the hip flexor synergy amplitude, τHE is the hip extensor timing onset, and FBKF is the amplitude of a synergy comprised of ST, ILf, and SA each activated similarly. A fixed amplitude flexor synergy of amplitude AKF = 0.30, with the muscle weights from the unadjusted simulation, was also combined in the simulation. As earlier, τHF = 0, τKF = 50 ms, AHE = 1.7. The muscle weights were given earlier.

We simulated wiping from 12 starting positions using this model. Note that the modulation of FBKF was exponential in this model and thus nonlinearly related to the polar coordinate distance. Figure 9 shows the starting positions relative to the hip joint (black dot) and the ankle paths from each starting position. The control structure produced both accurate (±4 mm from the central wiping target) and reasonably straight trajectories (mean path curvature = 1.3 ± 0.1 [SD]) across the workspace of the hindlimb. In addition, mean data for experimental frogs are shown when feedback is intact (Aff) and after deafferentation (Dff) in Table 1. These ranges closely match the variations of kinematic measures in spinal frogs in afferented conditions (Kargo and Giszter 2000b), as shown in Fig. 9B, and the fixed model matched the range of deafferented responses (Fig. 9B).

A: wiping simulations were performed at 12 different positions spanning the horizontal workspace of the hindlimb (marked 1–12). Lines represent the ankle paths of the model from each position, black dots represent the left and right hip joints, and the gray oval represents the target region where the wiping limb would normally contact the contralateral target limb. The control structure performed adequately in that ankle paths contacted the target region and were reasonably straight (and therefore direct). B: measures of trajectory errors (top) and path curvature (bottom) are compared graphically in 2 versions of the model and in intact and deafferented spinal frogs. The fixed motor pattern model matched the deafferented frog, whereas the fully adjusted model using the simple proprioceptive operations matched or exceeded the performance of the spinalized, but otherwise intact frog. See Table 1 for more detailed analysis of the model variants tested.

We next reexamined the composition of the hip flexor synergy. We examined ratios of muscles in synergies across our recordings of 24 frogs to see whether SA was especially variable in the hip flexor synergy. However, it was not variable, showing tight ratiometric couplings (see Fig. 10A). We therefore relaxed our initial minimum activation criterion for the hip flexor synergy, which had eliminated the SA muscle from participation in the initial simulations. We reasoned that the local isometric force description that we used in search and optimization might miss features of viscoelastic behavior, which were important in muscle choices used in real spinal synergies. The optimization of real spinal synergies might be more general in scope, extending beyond wiping alone and might represent more general dynamic and control features (e.g., see the 2D simulation and analysis of frog motion in Berniker et al. 2009). We restored SA with activation 0.5 to the HF synergy and rebalanced the synergy. Our analysis with this SA containing HF synergy showed that SA addition was as effective as or more effective than the previous HF synergy without SA (Fig. 10B). Thus as described previously (Giszter et al. 1993; Lemay et al. 2007; Loeb et al. 2000) similar isometric force pattern modularity could potentially be achieved with multiple synergies. The HF synergy with SA was very effective for proprioceptive adjustments across the workspace. Trajectories with path-ratio measures of straightness almost equivalent to Fig. 9 were obtained by altering only the HF amplitude and HE timing (i.e., using only Eqs. 13 and 14); there was thus less need to add the adjustments to KF that had been required in Eqs. 15 and 16 with the first version of HF synergy (e.g., Fig. 10, B and C shows trajectories without explicit KF regulation, average deviation 1.8, and average path length ratio 1.2). We also observed incidentally that the SA activation in the synergy appeared to enable wider speed ranges and brisker wipes, which were not part of the original simulation testing, and had not formed a performance criterion in selecting synergies. We observed that the HF synergy with SA added was sensitive to SA magnitude. In sensitivity analysis the highest deviated trajectories usually showed most SA variation (e.g., simulations in Fig. 10B), in keeping with the tight relationship with other muscles observed in the real frog (Fig. 10A).

Addition of SA to the HF synergy. A: to assess the stability of synergies across 24 frogs the ratios of muscle weights in the synergies extracted by infomax independent components analysis (Hart and Giszter 2004) were assessed for the 3 strongest muscles in each; muscles are stated in order above each plot. The ratios for 3 muscles form 2 paired ratios that define the third and this pair can be expressed as a point on a plane. We constructed these planes and the probabilities of observing points on the plane from the frogs' data. The thermogram “hot spots” or peaks represent high probability. Stereotyped synergies will have localized peaks. It can be seen SA occurs in highly localized peaks in each synergy in which it is present, including the HF synergy (labeled and highlighted). Synergies with RA as a member show most variability (middle panels). B: adding SA to the HF synergy with strength 0.5 and rebalancing the synergy gave similar simulated trajectories to the earlier simulation. HF was amplitude modulated with position and HE onset regulated with position; KF was fixed. The composition of the HF synergy was randomly varied by ≤50% of its muscle weights in the synergy to generate the spread of trajectories plotted. C: examples of effects of variation of all synergies by ≤50 and ≤100% of muscle weights in the synergies are shown from a sensitivity analysis; 5 sample trajectories are plotted starting from each location. D: coefficient of variation for the trajectory deviation (or distance to target) for muscle weight variations <30, 60, and 90% of the synergy weight. For <30% muscle variation coefficient of variation is <6%. It then rapidly rises to about 30% between 30 and 60% muscle weight variation. E: the path-ratio measure of straightness progressively increases, but remains low as synergies are altered. Path-straightness and ratio are dominated by initial configuration and the simple control regulations.

We then performed a more detailed sensitivity analysis to the HF and other synergies' muscle composition for wipes across the workspace. In the simulation sensitivity analysis all wipes had HF amplitude and HE timing adjusted using the simple spindle control mechanisms. This allowed us to explore the extent to which synergy muscle variations in our model degraded performance. We performed a sensitivity analysis in which individual muscle weights in their synergies could be randomly varied by any amount, up to a percentage ceiling level (with ceilings set at 10, 20, 30, 40, 50, and 100%). Figure 10C shows sample trajectories for muscle variations <50% and <100%. We then examined the variance and possible degradation of the performance measures (path ratio and deviation or distance to target/ overshoot) used earlier in the face of these synergy variations. For variations held to <30% of muscle activation, the trajectory behavior was fairly robust and the coefficient of variation (CV) of both these measures was <6% (Fig. 10, D and E). Trajectory parameters varied across the workspace much more as a result of the control of amplitude and timing of the synergies than as a result of the muscle composition variations in this range. CV values significantly and progressively increased as allowed variations in muscle balances were increased to reach peaks of roughly 30% for target deviations. Taken together, the sensitivity analysis showed that synergy muscle composition (e.g., presence of SA) was important for simpler control. However, the trajectories generated in this way with the synergies were robust for variations in the synergy composition <30% of the muscle magnitudes. The data generated by the simulation are thus compatible with robust control using simple algorithms achieved with pulsed synchronous synergies or primitives.

In this study, we tested the competence of modular adjustments and simple proprioceptive regulations within the framework of pulsed fixed synergies or primitives that has been identified in spinal frogs. We built a model of the frog limb and its wiping behavior that consisted of a kinematically realistic skeletomotor component (Kargo and Rome 2002; Kargo et al. 2002) combined with the primitive framework. We tested the individual, and the combined, capacity of the reflex adjustments in Kargo and Giszter (2008) to generate a range of adjusted trajectories using the three pulses seen in the placing trajectory of the wipe. Other primitives are engaged in later parts of wiping and in broader ranges of behavior (Hart and Giszter 2004; Kargo and Giszter 2000b), but were not necessary for the placing trajectory. The neural component included simple adjustments that regulated the amplitude and onset times of primitives forming the wiping motor pattern. Our study thus represents a test of the competence of a minimalist account of trajectory generation based on physiologically identified mechanisms.

There are several caveats with our study. Because of the redundancy of the muscular system, many synergies will produce the same endpoint force at a single limb position. Thus we might potentially have initially used incorrectly weighted synergies to model wiping. Indeed, the initial choice of a hip flexor (HF) synergy without the sartorius (SA), derived from the Monte Carlo and activation minimization procedure, was such a synergy. However, the ultimate muscle groups qualitatively matched synergies in Hart and Giszter (2004) and Kargo and Giszter (2000a, 2008). Our modeling across the workspace with the HF synergy without SA forced a revised composition of the controlled knee flexor group of muscles initially chosen. However, with SA the importance of this additional mechanism altering KF was significantly reduced. In a similar way, perhaps equifinal patterns might exist and have not been discovered by our Monte Carlo process. We did not seek to find patterns that demonstrated fully equifinal behavior, reaching the wiping target across the workspace without adjustment. However, such a pattern would be nonphysiological. After deafferentation it is known that the spinal frog overshoots the target region from many starting postures (Kargo and Giszter 2000b). In the simulations here, we tuned the amplitudes and onset times of synergies to minimize target misses and trajectory curvature. We showed that these tunings were readily obtained from simple modular pattern adjustments based on early feedback, similar to those described in Kargo and Giszter (2008). The final synergies represented a basis that was relatively insensitive to small variations and simulated trajectories were robust to changes in synergy balances when muscles were randomly deviated by <30% from their original values. However, intact frogs might optimize their movements based on many other criteria including energy expenditure, isochronous movement, and force minimization (e.g., see Kargo 2002). Intact frogs probably use feedback in several ways, extending beyond those simple controls presented here—e.g., they likely use feedback continuously during movement and can likely integrate the weak resting feedback. Further, muscle tone in spinal frogs is lower than that in intact frogs and accordingly spindle firing is likely higher in intact frogs even if quiescent. Clearly, it is almost certain that feedback from other populations of muscle spindles in addition to ST and ILf is used. Spindle populations may encode more complex state information and their feedback may be integrated in more complex ways than the simple arithmetic approaches used here (see Kargo and Giszter 2008; Weber 2007). This could aid in the exponential adjustment of knee flexion in our model. In addition, however, even in spinal frogs more sophisticated ways to estimate limb state beyond position are likely used, as they are in mammals (e.g., Scott and Loeb 1994; Weber et al. 2007).

Accurate kinematic models may be most crucial to the modeling we used (see Valero-Cuevas et al. 2003). We modeled only aponeurosis and wrapping structures as permitted in the SIMM environment, with possible drawbacks in some of the complex mechanisms of the frog leg. The adjustment algorithms tested here were examined with the Hill-type muscle models. Results might conceivably differ if more sophisticated muscle models were used (e.g., Cheng et al. 2000). However, we believe such differences will probably be quantitative rather than qualitative.

Although the controller used here was simple and omitted continuous on-line feedback effects after movement onset, the model captured the setup and execution of hindlimb wiping very well. Our data show that quite simple neural controls, which exploit motor primitives and reduce the degrees of freedom of the musculoskeletal plant, can perform well at adapting trajectories to variable external conditions, such as starting limb position. Our dynamic simulations based on primitives and simple controls captured the observed performance of real frogs. Similar primitive-like schemes have performed well at simulating human movements (Neptune et al. 2000; Raasch and Zajac 1999).

Our models, which more generally apply to spindle control effects in mammalian spinal cord, should probably be interpreted in light of the relations between frog and mammalian muscle spindles. Mammalian spindles provide information on muscle velocity and length in separate Ia and II afferent fibers and intrafusal muscle fibers within the spindle capsule can be independently controlled by gamma innervation and thus decoupled from activation of the muscle proper (Prochazka 1996). About 30% of mammalian spindles are also innervated by motoneurons controlling extrafusal muscle fibers (beta innervation), which directly couple intrafusal and extrafusal fibers (Illert et al. 1996). This beta innervation potentially supports positive feedback loops in which beta-induced increases in spindle firing act to increase the levels of beta activation (Houk 1972). In the frog, a single afferent carries both muscle velocity and length information and intrafusal innervation is purely beta (Eckhorn and Querfurth 1985). In frog spinal behaviors, our data suggest beta effects are substantial at EMG onset and before movement begins (e.g., as shown earlier in Murthy and Taylor 1978).

The early bursts of muscle spindle feedback allowed for accurate estimates of the starting limb position as the motor program began and could parameterize primitive-based algorithms for motor pattern construction. In the present study, it was difficult to assess whether variations in the premovement firing rates of ST and ILf spindles were better explained by variations in intrafusal fiber length or by variations in the level of beta drive. The frequency of the initial burst of muscle spindle firing to a constant-level burst of beta input depends on intrafusal muscle length (Corda et al. 1979). However, firing rates also depend on the level (and rate) of fusimotor activation (Matthews 1962). Furthermore, spindle feedback might act via a positive feedback loop to enhance beta levels and therefore further enhance spindle feedback (Houk 1972). Indeed we found that both the (mean) spike frequency of the premovement burst and the magnitude of the (rectified and integrated) homonymous EMG were linearly related to muscle length across positions. As a result of these interactions, ST and ILf spike frequency in this early period was related in a straightforward way to muscle length and, together, these spindles provided a coarse representation of the starting limb position. Polar estimates of ankle location were easily derived from arithmetic proprioceptive combinations of two muscle spindle sets. The use of addition and subtraction of firing rates of these two muscles' beta-spindle patterns to generate a limb representation could allow for simple low-level circuit mechanisms to readily parameterize the wiping behavior. Simplified coordinates of this type in locomotion and in spinal representation are suggested by the accumulated research of Bosco, Poppele, Borghese, Lacquaniti, Ivanenko, and colleagues (Borghese et al. 1996; Bosco and Poppele 2000, 2001; Bosco et al. 2006; Ivanenko et al. 2007).

Of course the simple control algorithms drawn from these representations could take several forms and need not directly use such a representation; moreover, much more information may be available in proprioception (see Weber et al. 2007). The modular adjustments observed in Kargo and Giszter (2008) that we tested in trajectory generation here might be considered an implicit low-level internal model of the behavior. The adjustment mechanisms did not explicitly represent the dynamic equations of motion or the transformational steps between a desired limb trajectory and muscle activation. However, they did successfully adjust the motor pattern in a simple fashion, to produce straighter and more accurate trajectories [i.e., to account for and therefore anticipate (external) mechanical variations at the different starting positions, e.g., variations in the force–length properties of muscles, moment arms, limb Jacobian, and variable interaction torques when moving from different starting positions]. There is biological value in a limb design and modular control that enables the use of simple algorithms and direct computations with sensors and that can implicitly anticipate and account for properties of the skeletomotor apparatus. Simpler and thus faster circuit networks can be used to implement them. Such simple modular control heuristics may be a principle of motor organization (McKay and Ting 2007). Using a simplified planar version of the model presented here, with fixed moment arms Berniker and colleagues (2009) were able to explore control strategies for a froglike limb in a more principled way. They showed that in the two-link planar case synergies like those observed in the frog and used here in the final simulations can capture more general features of the limb mechanics (their “natural dynamics synergies”) and thus can provide an optimal substrate for broader movement repertoires. In this study we tested the competence of the frog synergies and a set of physiologically identified simple controls in a 5 DOF simulation including all currently available anatomy. Shallow controls like those explored herein may also aid in the minor adjustments required for future learning—e.g., with growth or injury. However, if very simple use of proprioception such as that postulated here are indeed implemented in frogs biologically, we will need to identify their CNS cellular and network mechanisms (e.g., see Berkowitz 2001; Bosco and Poppele 2001; Prut and Fetz 1999).

The final model result and structure using the SA containing HF synergy remained consistent with pulsed primitives and fixed synergies and with our data from muscle vibration (Kargo and Giszter 2008) and from statistical decomposition (Hart and Giszter 2004). It is also worth noting that some of the changes needed across the workspace could in principle be implemented at the physiological level of the motoneuron (e.g., Hyngstrom et al. 2007).

The use of proprioception in the model here was postulated and modeled to occur briefly at movement onset. For rapid movements like wiping (<250 ms), continuous feedback control might not be appropriate because of relatively long transmission and electromechanical delays and the intermittency of the beta-driven feedback described here, which could potentially lead to instability or untimely corrections. An alternative to continuous feedback control is to preprogram muscle activation patterns that anticipate mechanical properties of the limb and environment for the upcoming movement (see Wolpert and Ghahramahni 2000). This strategy requires adequate information (available or predicted) about the current state of the limb and forces that will be experienced in the upcoming movement (e.g., Shadmehr and Mussa-Ivaldi 1994). In the present study, we provide preliminary evidence that spinal systems and reflex pathways might contribute implicitly to estimating limb position and setting up segmental motor patterns that minimize continuous feedback control (also see Bosco and Poppele 2000, 2001; Fukson et al. 1980). This form of control does not preclude the operation of other feedback controls that continuously regulate muscle and limb properties (Bonasera and Nichols 1996; Lin and Rymer 1998; Nichols and Houk 1976), EMG magnitude (Hiebert and Pearson 1999; Hyngstrom et al. 2007; McCrea 1998), and EMG timing (Grillner and Rossignol 1978; Hiebert et al. 1996). These might be especially important in perturbed movements, during fatigue, or immediately after musculoskeletal injury in animals with gamma motoneurons (e.g., see Whelan and Pearson 1997).

Several recent studies suggest that spinal networks can implement various computational algorithms (e.g., see Bosco and Poppele 2001; Fetz et al. 2000; Fukson et al. 1980; Pearson 2000; Prochazka 1996; Prut and Fetz 1999). We show here that rapid and simple trajectory adjustments can potentially be organized even at movement onset using fixed synergies activated in fixed duration pulses. These features along with module reusability will allow for the storage of larger behavioral repertoires in a finite neural space and may aid in bootstrapping a preliminary repertoire of motor behaviors.

In conclusion, our results show that very simple proprioceptive control algorithms—combined with modular motor pattern formation systems—can generate well-adjusted trajectories that match reflex kinematics organized in frog spinal cord and also capture many of the features of voluntary movements. These results provide some insight into how simple segmental control structures might contribute to reflex trajectory formation and more sophisticated voluntary limb movement control.

GRANTS

This work was supported by National Institute of Neurological Disorders and Stroke Grants NS-34640 and NS-40412 and National Science Foundation Grant IIS-0827684 to S. F. Giszter; National Institute of Arthritis and Musculoskeletal and Skin Diseases Grants AR-46125 and AR-38404 to L. C. Rome; National Institute of Child Health and Human Development Predoctoral Training Grant HD-07467 to Dr. Marion Murray; and an Allegheny Singer Research Institute graduate student grant to W. J. Kargo.

ACKNOWLEDGMENTS

We thank Dr. Rahman Davoodi for assistance with some MSMS use and Matlab implementations and Dr. Michel Lemay, J. Scabich, and the reviewers for helpful comments.

.
Towards a realistic biomechanical model of the thumb: the choice of kinematic description may be more critical than the solution method or the variability/uncertainty of musculoskeletal parameters. J Biomech36:
1019–1030, 2003.