- Email: [email protected]
Roles of the elements of the triphasic control signal
EXPERIMENTAL
NEUROLOGY
90,619-634
(1985)
Roles of the Elements of the Triphasic Control Signal BLAKEHANNAFORDANDLAWRENCESTARK' Department
of Electrical University Received
April
Engineering
of California.
and Computer Science, 483 Minor Berkeley* Calfornia 94720
5, 1985; revision
received
August
Hall,
5. 1985
In fast (time-optimal) movements about many joint systems, the triphasic EMG pattern has been observed. Although the first agonist burst obviously initiates the movement, the roles of the second and third bursts, appearing in the antagonist and agonist respectively, have been less clear. In this study, the timing of experimentally measured EMG signals led to construction of a three-pulse control signal that produced an accurate simulation of experimentally measured time-optimal head rotations using a sixth-order nonlinear model in conjunction with an optimization algorithm. By ablating pulses from the model control signal and observing the resulting dynamics, the roles of the three pulses can be assessed.As a result, the pulses can be designated PA, the action pulse (for the first agonist burst), PB, the braking pulse (for the antagonist burst), and PC, the clamping pulse (for the second agonist burst). Comparison of dynamic parameters from the simulated movements revealed strategies used to generate control Signak for movements of various speeds. 0 1985 Academic press, IX.
INTRODUCTION Head movements have many modes of control, many types of interaction with other movement systems, and many different trajectory types. Among the types that have been studied are fast, time-optimal, neurologically ballistic horizontal rotations. Time optimal head movements are those in which a motivated human subject is instructed and attempts to change the position of his head “as fast as possible” from one position to another. Neurologically ballistic movements (often incorrectly termed simply “ballistic”) are movements that are launched toward a known target with a predetermined initial Abbreviations: exaj, ex,,-agonist, antagonist excitation: kge-equivalent kilograms: PA, PB, PC-action, braking, clamping pulse. ’ We are pleased to acknowledge support from the National Institutes of Health’s Training Grant in Systems and Integrative Biology, from the Jet Propulsion Laboratory of the California Institute of Technology, and from the Essilor/Multioptics Corp., Paris and CA. 619 0014-4886/85 $3.00 Copyright Q 1985 by Academic Press, Inc. All rights of reproduction in any form reserved.
620
HANNAFORD
AND STARK
impetus. Unlike a truly ballistic trajectory, such as that of a thrown object, muscle forces are important throughout the movement. Although evoked in the laboratory under artificial conditions, these movements are similar to those naturally occurring, fast head movements required to orient the sensory apparatus toward an event detected far from the midline. The position, velocity, and acceleration time functions of these movements have been measured and characterized by such dynamic parameters as peak velocity, and positive and negative peaks of acceleration (27). These properties of time optimal head rotations have been plotted in the Main Sequence diagram [a log-log plot of dynamic parameters vs. magnitude; see (27)] to show an increasing relationship between magnitude and movement duration, peak velocity, and peak acceleration. Control Signals and Electromyograms. In analyzing movements of this type, records of the electromyogram (EMG) are commonly used to indicate timing and to some extent magnitude (12, 19) of activation of muscles involved in the movement. When fast movements, made in response to step target displacements, are analyzed, the EMG pattern recorded from antagonist muscles shows three pulses: an initial agonist burst, a burst of activity in the antagonist, and a final agonist burst. This triphasic EMG pattern has been observed in fast flexionextension movements of the elbow (2, 24, 25), in rotation of the wrist (17, 2 I), and in horizontal rotations of the head (10, 30). Ghez and Martin (7) observed it in fast leg movements in the cat. We denote these three EMG pulses PA, PB, and PC for the first agonist burst, the first antagonist burst, and the second agonist burst, respectively. Initially, the A, B, C notation will simply denote their temporal sequence. From the results of simulations, we will then develop mnemonic definitions. To some extent, the roles of the three pulses of activation are intuitively obvious. For example, it is clear that the initial agonist pulse which we call PA is responsible for initiating the movement and that the antagonist burst PB stops the movement, appearing when the movement is performed too fast for breaking by passive forces (16). The role of the second agonist burst PC has been less clear. Ghez and Martin, (7) speculated that it may have a viscous damping function, dissipating any excess energy remaining in the system after attainment of final target position. This paper is an attempt to precisely define the role of each of the three components of the classical triphasic response. This is achieved through simulation of a computer model of the muscles and plant driven by control signals based on the timing of recorded EMG activity. Modeling. Zangemeister et al. (28, 29) simulated the muscles and plant of the head movement system with a sixth-order nonlinear model incorporating Hill’s force-velocity relationship (1 l), two antagonist muscles, and an overdamped second-order plant (Fig. 1). Their model roughly matched experi-
TRIPHASIC
CONTROL
621
OF HEAD MOVEMENT
L \
fi J
\ \
K SL
/ KSR
\ \ \ \ \
‘+$vR li-
7-l KPL
KPR
\ LEFT
)
/ / / / / / / / /
RIGHT
FIG. 1. Sixth-order nonlinear model. Used to simulate experimentally measured time-optimal horizontal head rotations, the model consists of two identical antagonistic muscle models driving a second-order plant.
mentally measured Main Sequence dynamic peaks when driven by heuristically derived rectangular control signals. Versions of this model have had a fruitful history of application to many different physiologic systems. Stark (22) proposed and Atwood et al. (3) simulated a two-muscle model for understanding neurological control mechanisms. Cook and Stark (6), and Clark and Stark (5) used more detailed versions with appropriate parameter values to model saccadic and other eye movements, and it has been applied to the eyelid in modeling the dynamics of the blink ( 13). Optimization has been used in these modeling studies (26, 28) to estimate parameter or control signal values by solving for those values that minimize some performance criterion. In this study, an optimization algorithm was used to determine the amplitudes of neurologic control signal pulses needed to drive the model to match actual movements. In this case, the performance criterion was the root mean square error between model output and experimentally recorded movement dynamics. By minimizing this quantity, the optimization algorithm could find those pulse heights that result in a best fit to experimentally recorded data. METHODS
Experimental. Horizontal head position was measured by a precision potentiometer attached to a bicycle helmet frame worn by the subject. To allow for vertical head movement and subjects of differing heights, the potentiometer was coupled to the head through an adjustable fitting and universal joint.
622
HANNAFORD
AND STARK
Electromyographic activity was measured differentially with two Con-Med self-adhesive Ag-AgC1 disposable electrodes placed approximately 5 cm apart on the skin along the major axes of the left and right splenius capitus muscles. EMG and position data were digitized at 1000 Hz by an LSI 11-23 computer with 12-bit analog to digital converters (10). The subject’s head movements were made in response to light emitting diode (LED) targets alternately flashing at points on a curved screen 1 m from the subject’s head. The subject was made aware of the exact position of his head by a small spot of light projected from his helmet onto the screen. When the target illumination alternated between the two positions (at intervals of 4 s), the subject performed 20, 40, or 60” horizontal head movements. The subjects were instructed to move their heads “as fast as possible” to produce an intent to respond to the target in a time-optimal manner. Selection of the fastest five movements resulted in a stereotyped ensemble of movements which could be aligned in time (according to the zero crossing of acceleration) and averaged along with their rectified EMGs. More than seven subjects have been analyzed. The data presented are an ensemble average of five 40” movements from a single, typical subject. A more detailed discussion of ensemble averaging, selection, time alignment, etc. is presented in (10). Simulation results for 20 and 60” movements were also computed and agreed qualitatively with 40” results. Modeling and Optimization. For simulation of the horizontal head rotation system, we used the sixth-order nonlinear model developed by Zangemeister et al. (28, 29). This model consists of two identical, antagonistic muscle elements driving a second-order plant (Fig. 1). The muscle elements have a force generator driven through a time delay of 20 to 50 ms and a first-order, lowpass filter representing the calcium activation process. The control signals, exaB and exat (standing for agonist and antagonist excitation), ranged from zero to one to represent the possible range of excitation from none to full excitation. To help the reader’s intuition, this signal can also be described in equivalent kilograms (kge) to suggest the steady-state force that would result from constant excitation at a given level. In parallel with the force generator is a nonlinear viscous element representing A. V. Hill’s force-velocity relationship (11). Force was transmitted to the load through a series elastic element representing the properties of muscle tendon and attached cross bridges. This system is modeled and simulated by a set of six state equations and two ancillary equations (Table 1). The values used for the parameters (Table 2) are based on previous work (28, 29) and recent, improved estimates (J. M. Winters, private communication). Control Signals. Control signals for modeling represent levels of excitation delivered to muscle by the nerve. This is plotted as a dimensionless value from zero to one representing an amount of nervous excitation (a product of
TRIPHASIC
CONTROL
OF HEAD MOVEMENT
623
TABLE I Equations for the Sixth-Order Nonlinear Model of Horizontal Head Rotation State equations \’ zz ,’ v
= HTL
-
&XL
-
4
L BL
c, = -HTR R (I = -0
- Ks(sR - s) BR
x - Bp v + Ks(sL - s) + Ks(xR - x) J
dhl, = ~X.B - HTL TA d/.,1,. z exat
- HTR TA
Ancillary equations ‘-1.25 HTL bh + vL
VL > 0
BL = .,
1.25 HTL -
hh
~(!2$y R
VL < 0
VR < 0
BR =
1.25 HTR bh
VR > 0
firing rate and recruitment) ranging from none to a maximum. One, the maximum excitation level, represents a neurological command calling for maximal muscle output. The excitation control signal is the beginning of a chain of model elements and signals which partially describe the dynamics of muscle [a more detailed treatment is available in ( 14)]. After being filtered by the first-order model of the activation process, the excitation signal becomes a signal representing activation level. The activation level is converted to a hypothetical tension (HTL/R). This hypothetical tension, in the force generator, is two steps removed from the load. First, force resulting from the product of the nonlinear viscous element and the velocity of the muscle node ( vLIR) is subtracted (representing Hill’s force-velocity relationship). The remaining force is applied
624
HANNAFORD
AND STARK
TABLE 2 Model Parameter Definitions and Values Name Muscle model Activation time const Hill’s constant: b Series elasticity Second-order plant Rotational inertia Parallel elasticity Parallel viscosity Control signals Agonist excitation Antagonist excitation
Symbol
TA bh KS J KP BP
exss exat
Value 50.0 350.0 350.0
Units ms deg s-’ gr-f deg-’
0.18 2.0 2.0
gr-f deg-‘s’ gr-f deg-’ gr-f deg-‘s
-
dimensionless dimensionless
-
to the load through the series elastic element which represents the properties of muscle tendon and attached cross bridges. Thus, although it is intuitively convenient to think of the excitation control signal in units of force, with maximum excitation corresponding to a force level of say, 60 kge, it should be kept in mind that a control signal of 0.5 will result in the corresponding 30 kg of output force only in the steady state when all activation, muscle, and load dynamics have died out. Optimization. The optimization algorithm is based on the method of Bremmerman (4) and used extensively by others (15,23,26, 28). Full details are available in Bremmerman’s paper and especially (15). The basic idea is to take a two-dimensional section of the optimization space along a line of random orientation and compute five values of the quantity to be minimized (in our case RMS error between model and experiment). A polynomial is then fitted to the five points from which a projected minimum is computed. The function is then evaluated at this projected point which is often a new minimum. The procedure is repeated at different random orientations until no further improvement is obtained. RESULTS Fit of Control Signal to the Electromyogram. In order to use the model to investigate the role of the three observed pulses of EMG activity, a control signal must be obtained with which to drive the model. Unfortunately, this signal is not directly measurable. However, a reasonable approximation to this signal is a series of rectangular pulses with switching times derived from
TRIPHASIC
CONTROL
OF
HEAD
MOVEMENT
625
the pulses observed in the EMG signal. Switching times were determined by visually approximating the time that each pulse attained 50% of its peak value. Although there is reason to believe that the actual control signals are shaped differently (9) we felt that rectangular control pulses would yield sufficient accuracy to answer the question at hand. In addition, although the exact relationship between EMG and muscle force under dynamic conditions is not known, it is generally agreed that while the EMG activity is high, a large portion of the muscle is active, and vice versa. Thus, rectangular pulses (of as yet unspecified amplitude) whose timing is based on that of EMGs are the most conservative estimate of the underlying neurologic control signal. The pulse switching times used were determined by visual inspection of the EMG signals (Fig. 2). Finding Pulse Heights and the Values Found. To completely specify these approximate control signals, the pulse amplitudes as well as the pulse switching times must be determined. To do this, we used the method of optimization described above to determine the three amplitude values (for the three pulses, PA, PB, and PC) which minimized the RMS difference between the recorded dynamics and the model’s output. The optimization process varied the three pulse heights (the pulse switching times were fixed) and applied the pulses as inputs to the model, computing mean squared error between the resulting model output and the experimentally recorded movements. Optimal Puke Heights. In the first optimization experiment, the three pulse heights were allowed to vary. The resulting optimal pulse heights had PC saturated at its minimum value: a two-pulse signal. Optimal Three-Pulse Fits. Because of the EMG evidence for a third pulse, a second optimization experiment was conducted in which PC was permitted to vary only above a certain minimum value. Although the resulting fits to experimental data were suboptimal, the resulting dynamics were still a very good fit to the experimental data (Fig. 2). The controller signals resulting from this experiment had PC equal to its minimum permissable amplitude. Because three distinct EMG pulses were clear in the experimental data, this threepulse control signal was more consistent with the EMG evidence (Table 3). Simulation results; Pulse roles. When the complete, three-pulse control signal was used as model input, the output model dynamics closely matched an experimentally recorded 40” movement (Fig. 3, trace ABC, see also Fig. 2). The role of each pulse in the movement was revealed by simulations in which pulses were selectively removed from the control signal (Fig. 3, traces AB and A). Trace AB is the model output resulting from a control signal consisting of only pulses A and B. The AB trace was identical with the ABC trace until 2 15 ms after the onset of PA. The resultant position trace attained the desired final position, but then drifted away from target position. This
626
HANNAFORD
AND STARK
s
‘ic a75 :: w 1.00
t
FIG. 2. Head movement dynamics. Experimental records (solid lines), and simulation of same movement (dashed lines). Plotted are position (degrees), velocity (O/s-‘), and acceleration (O/s-*) for the experimentally measured (solid line) and simulated (dashed line) movement. Also plotted are EMGs from agonist (ag) and antagonist (ant) muscles (I mV full scale) and rectangular approximations to the EMG bursts used as inputs to the model. EMGs and rectangular control signals contained three pulses: action pulse (PA), thk first agonist burst, braking pulse (PB), the antagonist burst, and clamping pulse (PC), the second agonist burst. Experimental data represent ensemble average of five movements having the same dynamics.
TRIPHASIC
CONTROL
627
OF HEAD MOVEMENT
TABLE 3 Three-Pulse Control Signal Values which Resulted in Model Output that Best Matched Experimentally Recorded Head Movementsa Pulse height (excitation) Movement amplitude (deg)
A
B
C
20” 40” 60”
0.14 0.77 0.92
0.18 0.43 0.83
0.03 0.08 0.08
’ Pulse switching times were determined from EMGs recorded from the same movements. Pulse heights are expressed in units of activation with I .O corresponding to 60 keg; see text.
negative drift was also evident on examination of the velocity trace in which the negative velocity phase had been extended in duration in the absence of PC. Trace A is the result of a simulation in which only the PA pulse was supplied. The resulting trajectory was a movement of 64”, approximately 150% of the intended (target) amplitude. The peak velocity attained by this movement of 609 O/s was the same as that of the complete (ABC) movement, but the positive velocity phase was extended in the absence of PB, and the velocity remained positive throughout the extent of the movement. This corresponded to the absence of a negative phase of the acceleration curve in the absence of PB. These results suggested that the three control signal pulses can be labeled A, B, C to suggest the functions Action, the initial pulse that started the movement, Braking, the antagonist pulse that initiated deceleration, and Clamping, the final agonist pulse that fixed final position, Simulation results; Constant amplitude. A similar set of simulations could be computed in which the movement amplitude was held constant (Fig. 4). In these simulations, the pulse heights or excitation levels of the pulses had to be recomputed to preserve the correct final state (defined as position = final target position, and velocity = 0.0). Trace ABC is the same three-pulse simulation as in Fig. 3. It represents the best three-pulse fit to experimentally measured time-optimal movements. When restricted to two control signal pulses, the optimization process found an excitation level of 0.098 for PB when PA was constrained to be 0.5 (excitation levels were referred to 1.0 = 60 kge). The resulting position trace attained final position with a peak velocity of only 493”/s compared with 609”/s in the ABC case and showed no overshoot; i.e., the velocity returned to zero with little or no negative phase.
628
HANNAFORD
AND STARK
ABC
I 0.0
I 0.1
I 0.2
I 03
I 0.4
I 0.5
Time Irl
FIG. 3. Ablation of later control pukes without altering earlier pulses. Roles of the three control signal pulses are revealed in simulations in which pulses have been removed from control signal (scales same as Fig. 2). ABC-complete three-pulse control signal (same as that of Fig. 2; note excellent fit to experimental data). AB-simulated movement in which PC, clamping pulse, has been removed, (note inability to maintain final position after movement). A-simulation in which only PA is present (resulting movement is approximately twice intended magnitude).
When only the action pulse was present, a still lower excitation level (0.26) was required to attain correct final state (Fig. 4, a). The resulting position trajectory was a slowly rising and slowly converging path to final position. The velocity had a peak value of only 347”/s and a slow decay to zero with no negative extent. Note the qualitative similarity between this simulation and the PA simulation of a larger-amplitude movement (Fig. 3, A). Main Sequence. Key movement parameters from the simulated movements are abstracted in the Main Sequence diagram; a plot of extrema and timing of those parameters (Fig. 5). Points in each of the plots, representing the
TRIPHASIC
I 0.0
CONTROL
I 0.1
I 0.2
OF HEAD MOVEMENT
I 0.3
I 0.4
629
I 0.5
Time II)
FtG. 4. Ablation of control pulses while constraining final movement amplitude. Same simulation experiment as in Fig. 2, but control signal pulse amplitudes have been recomputed to force constant movement amplitude. ABC-same as in Figs. 2 and 3. ol@-optimal two-pulse fit. aoptimal single-pulse movement. Note decrease in time to reach target, increase in peak velocity, etc., as pulses were added to control signal.
previously described simulated movements, appeared in horizontal and vertical lines (for numerical values see Table 4). In the horizontal lines, the movement magnitude varied, but the movement dynamical parameters, i.e., peak velocity, were relatively constant as pulses were subtracted from the simulated control signal. In the vertical lines, movement magnitude was constant but movement parameters varied as pulses were subtracted and the heights of the remaining ones were varied. For one movement parameter, O)min (peak stopping acceleration), the addition of control signal pulses did not cause a purely horizontal or vertical set of Main Sequence points. In this case, PB caused large increases in the magnitude of Omin as well as a change in movement magnitude.
630
HANNAFORD
AND STARK
10.
A+
. :. .:.‘:I
A. : ;t., 0
c: id,
0
0
v
100
MAGNITUDE,
deg
FIG. 5. Main Sequence diagram for simulated movements. Compare with earlier measured results (dots) ( IO). Plotted are peak velocity (V), peak positive acceleration (A+), and peak stopping acceleration (A-). Ablation of control signal (as in Fig. 2) caused horizontal shiR on the main sequence diagram. When movement magnitude was held constant (as in Fig. 3) points lie in a vertical line. Because of the compression inherent in a log plot, vertical shifts have been exaggerated for clarity. For actual dynamic peaks of simulated movements, see Table 4.
DISCUSSION Previous modeling studies of these time-optimal, horizontal head rotations (28,29) have produced movements with the magnitudes, peak velocities, and accelerations characteristic of experimentally recorded movements, matching experimentally obtained Main Sequence parameters. However, the timing of their dynamic peaks failed to match experimental data. Similarly, the timing of their control signals was not based on EMG switching times so their control signal pulses do not correspond to the apparent muscle activation patterns. Our study has added two additional constraints to the modeling and optimization process that have helped to eliminate these difficulties. By constraining the control signal to have pulse switching times that correspond to EMG pulses and by optimizing the model fit with respect to mean squared error between model and experimental trajectories, a more accurate simulation was obtained that allowed study of the effects of the three components of the triphasic controller signal.
TRIPHASIC
CONTROL
631
OF HEAD MOVEMENT
TABLE 4 Dynamical Peaks for Simulated Movements with and without Ablation of Control Signal Pulses Peak acceleration: Control signal
Magnitude 6%)
Peak velocity (deg/s)
ABC” AB A 4 a
40 40 64 40 40
609 609 609 493 347
Starting (deg/s*) 1.08 x 1.08 x 1.08 x 8.76 x 6.24 x
IO4 IO4 IO4 IO3 IO’
Stopping (de&) -1.47 -1.47 -5.08 -9.14 -3.23
x x X x x
lo4 IO4 IO3 103 IO3
’ ABC indicates complete triphasic control signal; AB, ablation of second agonist pulse, PC; A, ablation of all but initial agonist burst.
“Horizontal” and “Vertical” Strategies. The main sequence diagram of simulated movements (Fig. 5) suggests two alternative strategies for generating control signals for movements of various magnitudes and trajectories. The main sequence of experimental head movement data shows a general trend toward greater speed as movement magnitude increases. In the “horizontal strategy” (constant peak velocity), a faster movement of a given magnitude is achieved by selecting an initial agonist pulse (PA) which starts a movement having a peak velocity characteristic of a single-pulse movement to a larger magnitude. PB is then selected to truncate the movement at the desired amplitude, but doing this (perhaps because of activation and deactivation lags) calls for an amount of stored energy which would result in an error in final position; this is corrected by PC. The dynamics resulting from the three control signals ABC, AB, and A lie in a horizontal line in each of the movement parameters shown in the main sequence. In the “vertical strategy” (constant amplitude), movement magnitude is held constant and pulse heights are varied in order to reach the target in different elapsed times. For the slowest movements, a single agonist pulse (model run a) is sufficient to move the limb to final position. Above a certain velocity, a correspondingly larger single agonist pulse would result in overshooting the target and so an antagonist pulse is added to the control signal (model run alpha-beta). Lestienne showed this experimentally in the elbow flexion-extension system (16), and recent results of Marsden’s group (18) demonstrate this type of strategy in finger movements [Fig. 3 of ( 1S)]. Finally, for time-optimal movements, employing maximal muscle forces, three pulses are required (intersection of horizontal and vertical strategies; model run ABC).
632
HANNAFORD
AND
STARK
Another way to look at the vertical strategy is as “reverse programming” of control signal pulses. In this idea, the later pulses can be thought of as enabling the earlier ones. Specifically, PC allows PB to be larger without disturbing the movement’s final position and velocity. This in turn allows PA to be larger which thereby increases movement peak velocity. Thus, addition of later pulses of activity allows a greater PA which increases peak velocity (for time optimality) while preserving desired final state. An interesting insight into the dynamics of fast movements is provided by the one dynamic parameter, O)min,which does not move in a purely horizontal or vertical direction with the addition of the braking pulse. The addition of PB to the control signal adds force which contributes to stopping the motion of the head. However, passive viscous forces in the plant and more importantly in the muscle, also contribute significant stopping force. Thus, the addition of PB decreases movement amplitude in addition to increasing peak stopping acceleration. This increase results from the addition of active to passive breaking forces. Because our subjects had practice, were well motivated, and only their fastest movements were selected for analysis, the movements we present are time-optimal at least to the extent that this can be defined in a biological system. Thus, comparison of dynamics resulting from alternate control signals with the optimal dynamics resulting from the ABC control signal reveals the effects on time optimality of control signal changes. The vertical, or reverse programming strategy, then, is a way to vary the speed of a given magnitude movement by adding control signal pulses. However, the degree of optimality is determined in advance of the movement, because the later pulses are determined by the size of PA and the final desired state. Development of the Previous Model. In conducting this study, the use of the RMS error criterion added an additional requirement of accuracy to the simulation process. To achieve accurate simulation of movement dynamics, delay elements between excitation control signal and hypothetical excitation input (not present in earlier versions of the model) were found to be necessary and indeed are suggested by experimental evidence (1, 8). Limitations of the Techniques.In choosing to approximate activation commands by rectangular pulses, we have taken a conservative approach to the interpretation of EMG envelopes. It may well be that the shape of the pulses contains additional information that allows more accurate reproduction of experimentally observed dynamics. Indeed, the rectangular control signals we derived from EMG envelopes were dependent to some extent on a threshold by which we determined the on and off times for each pulse. This source of subjectivity can be eliminated by using the EMG envelopes themselves as the basic shape of the control signals or by obtaining computed control signals by the inverse modeling technique (9). The use of EMG envelopes or even
CUMULATIVE region CA3, mossy fiber-induced excitation: long-term potentiation components (guinea pig), 90, 529 subiculum age-related dendritic changes branching pattern complexity loss (rhesus monkey), 87,4 I2 synaptic density decline (rhesus monkey), 87,403 dendritic changes after aluminum tartrate injections (rabbit), 89, 530 postnatal synaptogenesis after in titer0 halothane exposure (rat), 89, 520 d-tubocuramine direct injection, selective neurotoxicity (rat), 89, 172 Histamine storage in isolated superior cervical ganglia (mouse.), 90,36 Histochemistry cytochrome oxidase activity assessment for monitoring ischemic muscle injury (rat), 88,265 Histopathology morphometric analysis of spinal cord alterations after contusion (rat), 88, 135 Hoffman reflex facilitative effect of a pure tone (human), 89, 592 Hormones adrenocortical, and brain growth: reversibility and differential sensitivity during development (rat), 90,44 Horseradish peroxidase injection into hypoglossal nerve, peripheral course of lingual motoneuron axons (rat), 87, 20 intramuscular administration, peroxidase activity at nodes of Ranvier (mouse, rabbit, rat), 90, 272 retrograde axoplasmic transport applied to cut end of greater splanchnic nerve, sympathetic neuron projections (rabbit), 87, 334 lingual muscle innervation by hypoglossal and facial nerves (cat), 87, 578 Hurler’s syndrome cortical dendrite structure and development, Golgi and computer morphometric analysis, 88, 652 Hydrocephalus cerebrospinal fluid pressure and resistance to
SUBJECT INDEX
723
absorption during development (mouse), 98, 162 6-Hydroxydopamine effect on convulsive thresholds, comparison with eboracin (mouse), 89, 1 -induced locus ceruleus lesion, effect on cerebellar cortex development (rat), 88, I50 intracerebroventricular, effect on catecholamine& fiber regeneration (rat hypothalamus), 88, 336 neonatal administration, role of catecholamines in mediating acquisition of amygdala-kindled seizures (rat), 90,588 Hydroxyurea DNA synthesis inhibition in stretched patagialis muscle, effect on other growth characteristics (chicken), 87, 487 Hypercapnia systemic, 2-deoxyghtcose uptake in central nervous system (peripherally chemodenervated rat), 88, 673 Hypematremic dehydration elimination by intravenous glucose and fructose solutions, effect on seizures and water electrolyte content (rabbit), 87, 249 Hypertension acute, induced by Aramin injection, bloodbrain barrier opening without brain edema induction (rabbit), 87, 198 transient, phenylephrine-induced respiratory effects during sleep (cat), 90, 173 Hypoglossal nerve axons of lingual motoneurons, horseradish peroxidase study (rat), 87, 20 lingual muscle innervation, horseradish peroxidase study (cat), 87, 578 nucleus, evoked potentials and single unit response during photic stimulation of retina (rabbit), 98, 34 1 Hypoglycemia insulin-induced severe, with hypothermia, blood-brain barrier dysfunction (rat), 87, 129 Hypotension sciatic nerve blood flay autoregulation (rat), 90, 139 Hypothalamus arcuate nucleus, neuronal activity changes by electroacupuncture (rat), 87, 118 central catecholaminergic neuron regenera-
634
HANNAFORD
AND
STARK
7. GHEZ, C., AND J. H. MARTIN. 1982. Control ofrapid limb movement in the cat. III. Agonistantagonist coupling. Exp. Bruin Res. 45: 1 I S-125. 8. GORDON, J., AND C. GHEZ. 1984. EMG patterns in antagonist muscles during isometric contraction in man: relations to response dynamics. Exp. Brain Res. 55: 167- 17 I. 9. HANNAFORD, B., W. S. KIM, AND L. STARK. 1985. Inverse modeling reveals control signal in neurologically ballistic movements. Math. BioSci., in press. 10. HANNAFORD, B., M. H., NAM, V. LAKSHMINARAYANAN, AND L. STARK. 1984. EMG as controller signal with viscous load. J. Mofor Behav. 16: 225-274. 1I. HILL, A. V. 1938. The heat of shortening and dynamic constraints of muscle. Proc. R. Sot. Land.
126: 136-195.
12. HOGAN, N. 1976. A review of the methods of processing EMG for use as a proportional control signal. BioMed. Eng. 11: 8 l-86. 13. HUNG, G., F. HSU, AND L. STARK. 1978. Dynamics of the human eye blink. Am J. Optom. Physiol.
Opt. 54: 678-690.
14. LEHMAN, S. 1983. A Detailed Biophysical Model of Extraoccular Muscle. Doctoral dissertation, Department of Biophysics, Univ. of California, Berkeley. 15. LEHMAN, S. L., AND L. STARK. 1982. Three algorithms for interpreting models consisting of ordinary differential equations: sensitivity coefficientism, sensitivity functions, global optimization. Mud BioSci. 62: 107-122. 16. LESTIENNE, F. 1979. Effects of inertial load and velocity on the braking process of voluntary limb movements. Exp. Brain Res. 35: 407-4 18. 17. LITVINTSEV, A. I., AND N. S. SEROPYAN. 1977. Muscular control of movements with one degree of freedom. I. Single movements. Avfomnt. Telemekh. 5: 88-102. 18. MARSDEN, C. D., J. A. OBESO, AND J. C. ROTHWELL. 1983. The function of the agonist muscle during fast limb movement in man. Physiol. (London) 335: I-13. 19. PATLA, A. E., B. S. HUDGINS, P. A. PARKER, AND R. N. SCOTT. Myoelectric signal a quantitative measure of muscle mechanical output. Med. Biol. Eng. Comput. 20: 3 19-320. 20. PERRY,J., AND G. A. BEKEY. 1981. EMG-force relationships in skeletal muscle. CRC Crif. Rev. Biomed. Eng. 8: l-22. 2 I. SANES,J. N., AND V. A. JENNINGS.1984. Centrally programmed patterns of muscle activity in voluntary motor behavior in humans. Exp. Brain Res. 54: 23-32. 22. STARK, L. 1961. Neurological organization of the control system for movement. Q, Prog. Rep. No. 61, Res. Lab. Electronics. MIT, pp. 215-217, April 15. 23. SUN, F., W. KRENZ, AND L. STARK. 1983. A systems model for the pupil size effect. I. Transient data. Biol. Cybern. 48: 101-108. 24. WACHHOLDER, K., H. ALTENBURGER. 1926. Beitrage zur Physiologic der willkurlichen Bewegung. Pfltigers Arch. 214: 642-66 I. 25. WADMAN, W. J., J. J. D. VAN DER GON, AND R. J. A. DERKSEN. 1980. Muscle activation patterns for fast goal directed arm movements. J. Human Move. Stu. Vol. 6: 19-37. 26. WINTERS,J. W., M. H. NAM, AND L. STARK. 1984. Modeling dynamical Interaction between fast and slow movements: fast saccadic eye movement behavior in the presence of slower VOR. Math. BioSci. 68: 159. 27. ZANGEMEISTER, W. H., AND L. STARK. 1981. Dynamics of head movement trajectories: main sequence relationship. Exp. Neural. 71: 76-91. 28. ZANGEMEISTER,W. H., S. LEHMAN, AND L. STARK. I98 I. Sensitivity analysis and optimization for a head movement model. Biol. Cybern. 41: 33-45. 29. ZANGEMEISTER, W. H., S. LEHMAN, AND L. STARK. 1981. Simulation of head movement trajectories: model and fit to main sequence. Biol. Cybern. 41: 19-32. 30. ZANGEMEISTER,W. H., L. STARK, 0. MEIENBERG,AND T. WAITE. 1982. Neurological control of head rotations: electromyographic evidence. J. Neurol. Sci. 55: l-14.
NEUROLOGY
90,619-634
(1985)
Roles of the Elements of the Triphasic Control Signal BLAKEHANNAFORDANDLAWRENCESTARK' Department
of Electrical University Received
April
Engineering
of California.
and Computer Science, 483 Minor Berkeley* Calfornia 94720
5, 1985; revision
received
August
Hall,
5. 1985
In fast (time-optimal) movements about many joint systems, the triphasic EMG pattern has been observed. Although the first agonist burst obviously initiates the movement, the roles of the second and third bursts, appearing in the antagonist and agonist respectively, have been less clear. In this study, the timing of experimentally measured EMG signals led to construction of a three-pulse control signal that produced an accurate simulation of experimentally measured time-optimal head rotations using a sixth-order nonlinear model in conjunction with an optimization algorithm. By ablating pulses from the model control signal and observing the resulting dynamics, the roles of the three pulses can be assessed.As a result, the pulses can be designated PA, the action pulse (for the first agonist burst), PB, the braking pulse (for the antagonist burst), and PC, the clamping pulse (for the second agonist burst). Comparison of dynamic parameters from the simulated movements revealed strategies used to generate control Signak for movements of various speeds. 0 1985 Academic press, IX.
INTRODUCTION Head movements have many modes of control, many types of interaction with other movement systems, and many different trajectory types. Among the types that have been studied are fast, time-optimal, neurologically ballistic horizontal rotations. Time optimal head movements are those in which a motivated human subject is instructed and attempts to change the position of his head “as fast as possible” from one position to another. Neurologically ballistic movements (often incorrectly termed simply “ballistic”) are movements that are launched toward a known target with a predetermined initial Abbreviations: exaj, ex,,-agonist, antagonist excitation: kge-equivalent kilograms: PA, PB, PC-action, braking, clamping pulse. ’ We are pleased to acknowledge support from the National Institutes of Health’s Training Grant in Systems and Integrative Biology, from the Jet Propulsion Laboratory of the California Institute of Technology, and from the Essilor/Multioptics Corp., Paris and CA. 619 0014-4886/85 $3.00 Copyright Q 1985 by Academic Press, Inc. All rights of reproduction in any form reserved.
620
HANNAFORD
AND STARK
impetus. Unlike a truly ballistic trajectory, such as that of a thrown object, muscle forces are important throughout the movement. Although evoked in the laboratory under artificial conditions, these movements are similar to those naturally occurring, fast head movements required to orient the sensory apparatus toward an event detected far from the midline. The position, velocity, and acceleration time functions of these movements have been measured and characterized by such dynamic parameters as peak velocity, and positive and negative peaks of acceleration (27). These properties of time optimal head rotations have been plotted in the Main Sequence diagram [a log-log plot of dynamic parameters vs. magnitude; see (27)] to show an increasing relationship between magnitude and movement duration, peak velocity, and peak acceleration. Control Signals and Electromyograms. In analyzing movements of this type, records of the electromyogram (EMG) are commonly used to indicate timing and to some extent magnitude (12, 19) of activation of muscles involved in the movement. When fast movements, made in response to step target displacements, are analyzed, the EMG pattern recorded from antagonist muscles shows three pulses: an initial agonist burst, a burst of activity in the antagonist, and a final agonist burst. This triphasic EMG pattern has been observed in fast flexionextension movements of the elbow (2, 24, 25), in rotation of the wrist (17, 2 I), and in horizontal rotations of the head (10, 30). Ghez and Martin (7) observed it in fast leg movements in the cat. We denote these three EMG pulses PA, PB, and PC for the first agonist burst, the first antagonist burst, and the second agonist burst, respectively. Initially, the A, B, C notation will simply denote their temporal sequence. From the results of simulations, we will then develop mnemonic definitions. To some extent, the roles of the three pulses of activation are intuitively obvious. For example, it is clear that the initial agonist pulse which we call PA is responsible for initiating the movement and that the antagonist burst PB stops the movement, appearing when the movement is performed too fast for breaking by passive forces (16). The role of the second agonist burst PC has been less clear. Ghez and Martin, (7) speculated that it may have a viscous damping function, dissipating any excess energy remaining in the system after attainment of final target position. This paper is an attempt to precisely define the role of each of the three components of the classical triphasic response. This is achieved through simulation of a computer model of the muscles and plant driven by control signals based on the timing of recorded EMG activity. Modeling. Zangemeister et al. (28, 29) simulated the muscles and plant of the head movement system with a sixth-order nonlinear model incorporating Hill’s force-velocity relationship (1 l), two antagonist muscles, and an overdamped second-order plant (Fig. 1). Their model roughly matched experi-
TRIPHASIC
CONTROL
621
OF HEAD MOVEMENT
L \
fi J
\ \
K SL
/ KSR
\ \ \ \ \
‘+$vR li-
7-l KPL
KPR
\ LEFT
)
/ / / / / / / / /
RIGHT
FIG. 1. Sixth-order nonlinear model. Used to simulate experimentally measured time-optimal horizontal head rotations, the model consists of two identical antagonistic muscle models driving a second-order plant.
mentally measured Main Sequence dynamic peaks when driven by heuristically derived rectangular control signals. Versions of this model have had a fruitful history of application to many different physiologic systems. Stark (22) proposed and Atwood et al. (3) simulated a two-muscle model for understanding neurological control mechanisms. Cook and Stark (6), and Clark and Stark (5) used more detailed versions with appropriate parameter values to model saccadic and other eye movements, and it has been applied to the eyelid in modeling the dynamics of the blink ( 13). Optimization has been used in these modeling studies (26, 28) to estimate parameter or control signal values by solving for those values that minimize some performance criterion. In this study, an optimization algorithm was used to determine the amplitudes of neurologic control signal pulses needed to drive the model to match actual movements. In this case, the performance criterion was the root mean square error between model output and experimentally recorded movement dynamics. By minimizing this quantity, the optimization algorithm could find those pulse heights that result in a best fit to experimentally recorded data. METHODS
Experimental. Horizontal head position was measured by a precision potentiometer attached to a bicycle helmet frame worn by the subject. To allow for vertical head movement and subjects of differing heights, the potentiometer was coupled to the head through an adjustable fitting and universal joint.
622
HANNAFORD
AND STARK
Electromyographic activity was measured differentially with two Con-Med self-adhesive Ag-AgC1 disposable electrodes placed approximately 5 cm apart on the skin along the major axes of the left and right splenius capitus muscles. EMG and position data were digitized at 1000 Hz by an LSI 11-23 computer with 12-bit analog to digital converters (10). The subject’s head movements were made in response to light emitting diode (LED) targets alternately flashing at points on a curved screen 1 m from the subject’s head. The subject was made aware of the exact position of his head by a small spot of light projected from his helmet onto the screen. When the target illumination alternated between the two positions (at intervals of 4 s), the subject performed 20, 40, or 60” horizontal head movements. The subjects were instructed to move their heads “as fast as possible” to produce an intent to respond to the target in a time-optimal manner. Selection of the fastest five movements resulted in a stereotyped ensemble of movements which could be aligned in time (according to the zero crossing of acceleration) and averaged along with their rectified EMGs. More than seven subjects have been analyzed. The data presented are an ensemble average of five 40” movements from a single, typical subject. A more detailed discussion of ensemble averaging, selection, time alignment, etc. is presented in (10). Simulation results for 20 and 60” movements were also computed and agreed qualitatively with 40” results. Modeling and Optimization. For simulation of the horizontal head rotation system, we used the sixth-order nonlinear model developed by Zangemeister et al. (28, 29). This model consists of two identical, antagonistic muscle elements driving a second-order plant (Fig. 1). The muscle elements have a force generator driven through a time delay of 20 to 50 ms and a first-order, lowpass filter representing the calcium activation process. The control signals, exaB and exat (standing for agonist and antagonist excitation), ranged from zero to one to represent the possible range of excitation from none to full excitation. To help the reader’s intuition, this signal can also be described in equivalent kilograms (kge) to suggest the steady-state force that would result from constant excitation at a given level. In parallel with the force generator is a nonlinear viscous element representing A. V. Hill’s force-velocity relationship (11). Force was transmitted to the load through a series elastic element representing the properties of muscle tendon and attached cross bridges. This system is modeled and simulated by a set of six state equations and two ancillary equations (Table 1). The values used for the parameters (Table 2) are based on previous work (28, 29) and recent, improved estimates (J. M. Winters, private communication). Control Signals. Control signals for modeling represent levels of excitation delivered to muscle by the nerve. This is plotted as a dimensionless value from zero to one representing an amount of nervous excitation (a product of
TRIPHASIC
CONTROL
OF HEAD MOVEMENT
623
TABLE I Equations for the Sixth-Order Nonlinear Model of Horizontal Head Rotation State equations \’ zz ,’ v
= HTL
-
&XL
-
4
L BL
c, = -HTR R (I = -0
- Ks(sR - s) BR
x - Bp v + Ks(sL - s) + Ks(xR - x) J
dhl, = ~X.B - HTL TA d/.,1,. z exat
- HTR TA
Ancillary equations ‘-1.25 HTL bh + vL
VL > 0
BL = .,
1.25 HTL -
hh
~(!2$y R
VL < 0
VR < 0
BR =
1.25 HTR bh
VR > 0
firing rate and recruitment) ranging from none to a maximum. One, the maximum excitation level, represents a neurological command calling for maximal muscle output. The excitation control signal is the beginning of a chain of model elements and signals which partially describe the dynamics of muscle [a more detailed treatment is available in ( 14)]. After being filtered by the first-order model of the activation process, the excitation signal becomes a signal representing activation level. The activation level is converted to a hypothetical tension (HTL/R). This hypothetical tension, in the force generator, is two steps removed from the load. First, force resulting from the product of the nonlinear viscous element and the velocity of the muscle node ( vLIR) is subtracted (representing Hill’s force-velocity relationship). The remaining force is applied
624
HANNAFORD
AND STARK
TABLE 2 Model Parameter Definitions and Values Name Muscle model Activation time const Hill’s constant: b Series elasticity Second-order plant Rotational inertia Parallel elasticity Parallel viscosity Control signals Agonist excitation Antagonist excitation
Symbol
TA bh KS J KP BP
exss exat
Value 50.0 350.0 350.0
Units ms deg s-’ gr-f deg-’
0.18 2.0 2.0
gr-f deg-‘s’ gr-f deg-’ gr-f deg-‘s
-
dimensionless dimensionless
-
to the load through the series elastic element which represents the properties of muscle tendon and attached cross bridges. Thus, although it is intuitively convenient to think of the excitation control signal in units of force, with maximum excitation corresponding to a force level of say, 60 kge, it should be kept in mind that a control signal of 0.5 will result in the corresponding 30 kg of output force only in the steady state when all activation, muscle, and load dynamics have died out. Optimization. The optimization algorithm is based on the method of Bremmerman (4) and used extensively by others (15,23,26, 28). Full details are available in Bremmerman’s paper and especially (15). The basic idea is to take a two-dimensional section of the optimization space along a line of random orientation and compute five values of the quantity to be minimized (in our case RMS error between model and experiment). A polynomial is then fitted to the five points from which a projected minimum is computed. The function is then evaluated at this projected point which is often a new minimum. The procedure is repeated at different random orientations until no further improvement is obtained. RESULTS Fit of Control Signal to the Electromyogram. In order to use the model to investigate the role of the three observed pulses of EMG activity, a control signal must be obtained with which to drive the model. Unfortunately, this signal is not directly measurable. However, a reasonable approximation to this signal is a series of rectangular pulses with switching times derived from
TRIPHASIC
CONTROL
OF
HEAD
MOVEMENT
625
the pulses observed in the EMG signal. Switching times were determined by visually approximating the time that each pulse attained 50% of its peak value. Although there is reason to believe that the actual control signals are shaped differently (9) we felt that rectangular control pulses would yield sufficient accuracy to answer the question at hand. In addition, although the exact relationship between EMG and muscle force under dynamic conditions is not known, it is generally agreed that while the EMG activity is high, a large portion of the muscle is active, and vice versa. Thus, rectangular pulses (of as yet unspecified amplitude) whose timing is based on that of EMGs are the most conservative estimate of the underlying neurologic control signal. The pulse switching times used were determined by visual inspection of the EMG signals (Fig. 2). Finding Pulse Heights and the Values Found. To completely specify these approximate control signals, the pulse amplitudes as well as the pulse switching times must be determined. To do this, we used the method of optimization described above to determine the three amplitude values (for the three pulses, PA, PB, and PC) which minimized the RMS difference between the recorded dynamics and the model’s output. The optimization process varied the three pulse heights (the pulse switching times were fixed) and applied the pulses as inputs to the model, computing mean squared error between the resulting model output and the experimentally recorded movements. Optimal Puke Heights. In the first optimization experiment, the three pulse heights were allowed to vary. The resulting optimal pulse heights had PC saturated at its minimum value: a two-pulse signal. Optimal Three-Pulse Fits. Because of the EMG evidence for a third pulse, a second optimization experiment was conducted in which PC was permitted to vary only above a certain minimum value. Although the resulting fits to experimental data were suboptimal, the resulting dynamics were still a very good fit to the experimental data (Fig. 2). The controller signals resulting from this experiment had PC equal to its minimum permissable amplitude. Because three distinct EMG pulses were clear in the experimental data, this threepulse control signal was more consistent with the EMG evidence (Table 3). Simulation results; Pulse roles. When the complete, three-pulse control signal was used as model input, the output model dynamics closely matched an experimentally recorded 40” movement (Fig. 3, trace ABC, see also Fig. 2). The role of each pulse in the movement was revealed by simulations in which pulses were selectively removed from the control signal (Fig. 3, traces AB and A). Trace AB is the model output resulting from a control signal consisting of only pulses A and B. The AB trace was identical with the ABC trace until 2 15 ms after the onset of PA. The resultant position trace attained the desired final position, but then drifted away from target position. This
626
HANNAFORD
AND STARK
s
‘ic a75 :: w 1.00
t
FIG. 2. Head movement dynamics. Experimental records (solid lines), and simulation of same movement (dashed lines). Plotted are position (degrees), velocity (O/s-‘), and acceleration (O/s-*) for the experimentally measured (solid line) and simulated (dashed line) movement. Also plotted are EMGs from agonist (ag) and antagonist (ant) muscles (I mV full scale) and rectangular approximations to the EMG bursts used as inputs to the model. EMGs and rectangular control signals contained three pulses: action pulse (PA), thk first agonist burst, braking pulse (PB), the antagonist burst, and clamping pulse (PC), the second agonist burst. Experimental data represent ensemble average of five movements having the same dynamics.
TRIPHASIC
CONTROL
627
OF HEAD MOVEMENT
TABLE 3 Three-Pulse Control Signal Values which Resulted in Model Output that Best Matched Experimentally Recorded Head Movementsa Pulse height (excitation) Movement amplitude (deg)
A
B
C
20” 40” 60”
0.14 0.77 0.92
0.18 0.43 0.83
0.03 0.08 0.08
’ Pulse switching times were determined from EMGs recorded from the same movements. Pulse heights are expressed in units of activation with I .O corresponding to 60 keg; see text.
negative drift was also evident on examination of the velocity trace in which the negative velocity phase had been extended in duration in the absence of PC. Trace A is the result of a simulation in which only the PA pulse was supplied. The resulting trajectory was a movement of 64”, approximately 150% of the intended (target) amplitude. The peak velocity attained by this movement of 609 O/s was the same as that of the complete (ABC) movement, but the positive velocity phase was extended in the absence of PB, and the velocity remained positive throughout the extent of the movement. This corresponded to the absence of a negative phase of the acceleration curve in the absence of PB. These results suggested that the three control signal pulses can be labeled A, B, C to suggest the functions Action, the initial pulse that started the movement, Braking, the antagonist pulse that initiated deceleration, and Clamping, the final agonist pulse that fixed final position, Simulation results; Constant amplitude. A similar set of simulations could be computed in which the movement amplitude was held constant (Fig. 4). In these simulations, the pulse heights or excitation levels of the pulses had to be recomputed to preserve the correct final state (defined as position = final target position, and velocity = 0.0). Trace ABC is the same three-pulse simulation as in Fig. 3. It represents the best three-pulse fit to experimentally measured time-optimal movements. When restricted to two control signal pulses, the optimization process found an excitation level of 0.098 for PB when PA was constrained to be 0.5 (excitation levels were referred to 1.0 = 60 kge). The resulting position trace attained final position with a peak velocity of only 493”/s compared with 609”/s in the ABC case and showed no overshoot; i.e., the velocity returned to zero with little or no negative phase.
628
HANNAFORD
AND STARK
ABC
I 0.0
I 0.1
I 0.2
I 03
I 0.4
I 0.5
Time Irl
FIG. 3. Ablation of later control pukes without altering earlier pulses. Roles of the three control signal pulses are revealed in simulations in which pulses have been removed from control signal (scales same as Fig. 2). ABC-complete three-pulse control signal (same as that of Fig. 2; note excellent fit to experimental data). AB-simulated movement in which PC, clamping pulse, has been removed, (note inability to maintain final position after movement). A-simulation in which only PA is present (resulting movement is approximately twice intended magnitude).
When only the action pulse was present, a still lower excitation level (0.26) was required to attain correct final state (Fig. 4, a). The resulting position trajectory was a slowly rising and slowly converging path to final position. The velocity had a peak value of only 347”/s and a slow decay to zero with no negative extent. Note the qualitative similarity between this simulation and the PA simulation of a larger-amplitude movement (Fig. 3, A). Main Sequence. Key movement parameters from the simulated movements are abstracted in the Main Sequence diagram; a plot of extrema and timing of those parameters (Fig. 5). Points in each of the plots, representing the
TRIPHASIC
I 0.0
CONTROL
I 0.1
I 0.2
OF HEAD MOVEMENT
I 0.3
I 0.4
629
I 0.5
Time II)
FtG. 4. Ablation of control pulses while constraining final movement amplitude. Same simulation experiment as in Fig. 2, but control signal pulse amplitudes have been recomputed to force constant movement amplitude. ABC-same as in Figs. 2 and 3. ol@-optimal two-pulse fit. aoptimal single-pulse movement. Note decrease in time to reach target, increase in peak velocity, etc., as pulses were added to control signal.
previously described simulated movements, appeared in horizontal and vertical lines (for numerical values see Table 4). In the horizontal lines, the movement magnitude varied, but the movement dynamical parameters, i.e., peak velocity, were relatively constant as pulses were subtracted from the simulated control signal. In the vertical lines, movement magnitude was constant but movement parameters varied as pulses were subtracted and the heights of the remaining ones were varied. For one movement parameter, O)min (peak stopping acceleration), the addition of control signal pulses did not cause a purely horizontal or vertical set of Main Sequence points. In this case, PB caused large increases in the magnitude of Omin as well as a change in movement magnitude.
630
HANNAFORD
AND STARK
10.
A+
. :. .:.‘:I
A. : ;t., 0
c: id,
0
0
v
100
MAGNITUDE,
deg
FIG. 5. Main Sequence diagram for simulated movements. Compare with earlier measured results (dots) ( IO). Plotted are peak velocity (V), peak positive acceleration (A+), and peak stopping acceleration (A-). Ablation of control signal (as in Fig. 2) caused horizontal shiR on the main sequence diagram. When movement magnitude was held constant (as in Fig. 3) points lie in a vertical line. Because of the compression inherent in a log plot, vertical shifts have been exaggerated for clarity. For actual dynamic peaks of simulated movements, see Table 4.
DISCUSSION Previous modeling studies of these time-optimal, horizontal head rotations (28,29) have produced movements with the magnitudes, peak velocities, and accelerations characteristic of experimentally recorded movements, matching experimentally obtained Main Sequence parameters. However, the timing of their dynamic peaks failed to match experimental data. Similarly, the timing of their control signals was not based on EMG switching times so their control signal pulses do not correspond to the apparent muscle activation patterns. Our study has added two additional constraints to the modeling and optimization process that have helped to eliminate these difficulties. By constraining the control signal to have pulse switching times that correspond to EMG pulses and by optimizing the model fit with respect to mean squared error between model and experimental trajectories, a more accurate simulation was obtained that allowed study of the effects of the three components of the triphasic controller signal.
TRIPHASIC
CONTROL
631
OF HEAD MOVEMENT
TABLE 4 Dynamical Peaks for Simulated Movements with and without Ablation of Control Signal Pulses Peak acceleration: Control signal
Magnitude 6%)
Peak velocity (deg/s)
ABC” AB A 4 a
40 40 64 40 40
609 609 609 493 347
Starting (deg/s*) 1.08 x 1.08 x 1.08 x 8.76 x 6.24 x
IO4 IO4 IO4 IO3 IO’
Stopping (de&) -1.47 -1.47 -5.08 -9.14 -3.23
x x X x x
lo4 IO4 IO3 103 IO3
’ ABC indicates complete triphasic control signal; AB, ablation of second agonist pulse, PC; A, ablation of all but initial agonist burst.
“Horizontal” and “Vertical” Strategies. The main sequence diagram of simulated movements (Fig. 5) suggests two alternative strategies for generating control signals for movements of various magnitudes and trajectories. The main sequence of experimental head movement data shows a general trend toward greater speed as movement magnitude increases. In the “horizontal strategy” (constant peak velocity), a faster movement of a given magnitude is achieved by selecting an initial agonist pulse (PA) which starts a movement having a peak velocity characteristic of a single-pulse movement to a larger magnitude. PB is then selected to truncate the movement at the desired amplitude, but doing this (perhaps because of activation and deactivation lags) calls for an amount of stored energy which would result in an error in final position; this is corrected by PC. The dynamics resulting from the three control signals ABC, AB, and A lie in a horizontal line in each of the movement parameters shown in the main sequence. In the “vertical strategy” (constant amplitude), movement magnitude is held constant and pulse heights are varied in order to reach the target in different elapsed times. For the slowest movements, a single agonist pulse (model run a) is sufficient to move the limb to final position. Above a certain velocity, a correspondingly larger single agonist pulse would result in overshooting the target and so an antagonist pulse is added to the control signal (model run alpha-beta). Lestienne showed this experimentally in the elbow flexion-extension system (16), and recent results of Marsden’s group (18) demonstrate this type of strategy in finger movements [Fig. 3 of ( 1S)]. Finally, for time-optimal movements, employing maximal muscle forces, three pulses are required (intersection of horizontal and vertical strategies; model run ABC).
632
HANNAFORD
AND
STARK
Another way to look at the vertical strategy is as “reverse programming” of control signal pulses. In this idea, the later pulses can be thought of as enabling the earlier ones. Specifically, PC allows PB to be larger without disturbing the movement’s final position and velocity. This in turn allows PA to be larger which thereby increases movement peak velocity. Thus, addition of later pulses of activity allows a greater PA which increases peak velocity (for time optimality) while preserving desired final state. An interesting insight into the dynamics of fast movements is provided by the one dynamic parameter, O)min,which does not move in a purely horizontal or vertical direction with the addition of the braking pulse. The addition of PB to the control signal adds force which contributes to stopping the motion of the head. However, passive viscous forces in the plant and more importantly in the muscle, also contribute significant stopping force. Thus, the addition of PB decreases movement amplitude in addition to increasing peak stopping acceleration. This increase results from the addition of active to passive breaking forces. Because our subjects had practice, were well motivated, and only their fastest movements were selected for analysis, the movements we present are time-optimal at least to the extent that this can be defined in a biological system. Thus, comparison of dynamics resulting from alternate control signals with the optimal dynamics resulting from the ABC control signal reveals the effects on time optimality of control signal changes. The vertical, or reverse programming strategy, then, is a way to vary the speed of a given magnitude movement by adding control signal pulses. However, the degree of optimality is determined in advance of the movement, because the later pulses are determined by the size of PA and the final desired state. Development of the Previous Model. In conducting this study, the use of the RMS error criterion added an additional requirement of accuracy to the simulation process. To achieve accurate simulation of movement dynamics, delay elements between excitation control signal and hypothetical excitation input (not present in earlier versions of the model) were found to be necessary and indeed are suggested by experimental evidence (1, 8). Limitations of the Techniques.In choosing to approximate activation commands by rectangular pulses, we have taken a conservative approach to the interpretation of EMG envelopes. It may well be that the shape of the pulses contains additional information that allows more accurate reproduction of experimentally observed dynamics. Indeed, the rectangular control signals we derived from EMG envelopes were dependent to some extent on a threshold by which we determined the on and off times for each pulse. This source of subjectivity can be eliminated by using the EMG envelopes themselves as the basic shape of the control signals or by obtaining computed control signals by the inverse modeling technique (9). The use of EMG envelopes or even
CUMULATIVE region CA3, mossy fiber-induced excitation: long-term potentiation components (guinea pig), 90, 529 subiculum age-related dendritic changes branching pattern complexity loss (rhesus monkey), 87,4 I2 synaptic density decline (rhesus monkey), 87,403 dendritic changes after aluminum tartrate injections (rabbit), 89, 530 postnatal synaptogenesis after in titer0 halothane exposure (rat), 89, 520 d-tubocuramine direct injection, selective neurotoxicity (rat), 89, 172 Histamine storage in isolated superior cervical ganglia (mouse.), 90,36 Histochemistry cytochrome oxidase activity assessment for monitoring ischemic muscle injury (rat), 88,265 Histopathology morphometric analysis of spinal cord alterations after contusion (rat), 88, 135 Hoffman reflex facilitative effect of a pure tone (human), 89, 592 Hormones adrenocortical, and brain growth: reversibility and differential sensitivity during development (rat), 90,44 Horseradish peroxidase injection into hypoglossal nerve, peripheral course of lingual motoneuron axons (rat), 87, 20 intramuscular administration, peroxidase activity at nodes of Ranvier (mouse, rabbit, rat), 90, 272 retrograde axoplasmic transport applied to cut end of greater splanchnic nerve, sympathetic neuron projections (rabbit), 87, 334 lingual muscle innervation by hypoglossal and facial nerves (cat), 87, 578 Hurler’s syndrome cortical dendrite structure and development, Golgi and computer morphometric analysis, 88, 652 Hydrocephalus cerebrospinal fluid pressure and resistance to
SUBJECT INDEX
723
absorption during development (mouse), 98, 162 6-Hydroxydopamine effect on convulsive thresholds, comparison with eboracin (mouse), 89, 1 -induced locus ceruleus lesion, effect on cerebellar cortex development (rat), 88, I50 intracerebroventricular, effect on catecholamine& fiber regeneration (rat hypothalamus), 88, 336 neonatal administration, role of catecholamines in mediating acquisition of amygdala-kindled seizures (rat), 90,588 Hydroxyurea DNA synthesis inhibition in stretched patagialis muscle, effect on other growth characteristics (chicken), 87, 487 Hypercapnia systemic, 2-deoxyghtcose uptake in central nervous system (peripherally chemodenervated rat), 88, 673 Hypematremic dehydration elimination by intravenous glucose and fructose solutions, effect on seizures and water electrolyte content (rabbit), 87, 249 Hypertension acute, induced by Aramin injection, bloodbrain barrier opening without brain edema induction (rabbit), 87, 198 transient, phenylephrine-induced respiratory effects during sleep (cat), 90, 173 Hypoglossal nerve axons of lingual motoneurons, horseradish peroxidase study (rat), 87, 20 lingual muscle innervation, horseradish peroxidase study (cat), 87, 578 nucleus, evoked potentials and single unit response during photic stimulation of retina (rabbit), 98, 34 1 Hypoglycemia insulin-induced severe, with hypothermia, blood-brain barrier dysfunction (rat), 87, 129 Hypotension sciatic nerve blood flay autoregulation (rat), 90, 139 Hypothalamus arcuate nucleus, neuronal activity changes by electroacupuncture (rat), 87, 118 central catecholaminergic neuron regenera-
634
HANNAFORD
AND
STARK
7. GHEZ, C., AND J. H. MARTIN. 1982. Control ofrapid limb movement in the cat. III. Agonistantagonist coupling. Exp. Bruin Res. 45: 1 I S-125. 8. GORDON, J., AND C. GHEZ. 1984. EMG patterns in antagonist muscles during isometric contraction in man: relations to response dynamics. Exp. Brain Res. 55: 167- 17 I. 9. HANNAFORD, B., W. S. KIM, AND L. STARK. 1985. Inverse modeling reveals control signal in neurologically ballistic movements. Math. BioSci., in press. 10. HANNAFORD, B., M. H., NAM, V. LAKSHMINARAYANAN, AND L. STARK. 1984. EMG as controller signal with viscous load. J. Mofor Behav. 16: 225-274. 1I. HILL, A. V. 1938. The heat of shortening and dynamic constraints of muscle. Proc. R. Sot. Land.
126: 136-195.
12. HOGAN, N. 1976. A review of the methods of processing EMG for use as a proportional control signal. BioMed. Eng. 11: 8 l-86. 13. HUNG, G., F. HSU, AND L. STARK. 1978. Dynamics of the human eye blink. Am J. Optom. Physiol.
Opt. 54: 678-690.
14. LEHMAN, S. 1983. A Detailed Biophysical Model of Extraoccular Muscle. Doctoral dissertation, Department of Biophysics, Univ. of California, Berkeley. 15. LEHMAN, S. L., AND L. STARK. 1982. Three algorithms for interpreting models consisting of ordinary differential equations: sensitivity coefficientism, sensitivity functions, global optimization. Mud BioSci. 62: 107-122. 16. LESTIENNE, F. 1979. Effects of inertial load and velocity on the braking process of voluntary limb movements. Exp. Brain Res. 35: 407-4 18. 17. LITVINTSEV, A. I., AND N. S. SEROPYAN. 1977. Muscular control of movements with one degree of freedom. I. Single movements. Avfomnt. Telemekh. 5: 88-102. 18. MARSDEN, C. D., J. A. OBESO, AND J. C. ROTHWELL. 1983. The function of the agonist muscle during fast limb movement in man. Physiol. (London) 335: I-13. 19. PATLA, A. E., B. S. HUDGINS, P. A. PARKER, AND R. N. SCOTT. Myoelectric signal a quantitative measure of muscle mechanical output. Med. Biol. Eng. Comput. 20: 3 19-320. 20. PERRY,J., AND G. A. BEKEY. 1981. EMG-force relationships in skeletal muscle. CRC Crif. Rev. Biomed. Eng. 8: l-22. 2 I. SANES,J. N., AND V. A. JENNINGS.1984. Centrally programmed patterns of muscle activity in voluntary motor behavior in humans. Exp. Brain Res. 54: 23-32. 22. STARK, L. 1961. Neurological organization of the control system for movement. Q, Prog. Rep. No. 61, Res. Lab. Electronics. MIT, pp. 215-217, April 15. 23. SUN, F., W. KRENZ, AND L. STARK. 1983. A systems model for the pupil size effect. I. Transient data. Biol. Cybern. 48: 101-108. 24. WACHHOLDER, K., H. ALTENBURGER. 1926. Beitrage zur Physiologic der willkurlichen Bewegung. Pfltigers Arch. 214: 642-66 I. 25. WADMAN, W. J., J. J. D. VAN DER GON, AND R. J. A. DERKSEN. 1980. Muscle activation patterns for fast goal directed arm movements. J. Human Move. Stu. Vol. 6: 19-37. 26. WINTERS,J. W., M. H. NAM, AND L. STARK. 1984. Modeling dynamical Interaction between fast and slow movements: fast saccadic eye movement behavior in the presence of slower VOR. Math. BioSci. 68: 159. 27. ZANGEMEISTER, W. H., AND L. STARK. 1981. Dynamics of head movement trajectories: main sequence relationship. Exp. Neural. 71: 76-91. 28. ZANGEMEISTER,W. H., S. LEHMAN, AND L. STARK. I98 I. Sensitivity analysis and optimization for a head movement model. Biol. Cybern. 41: 33-45. 29. ZANGEMEISTER, W. H., S. LEHMAN, AND L. STARK. 1981. Simulation of head movement trajectories: model and fit to main sequence. Biol. Cybern. 41: 19-32. 30. ZANGEMEISTER,W. H., L. STARK, 0. MEIENBERG,AND T. WAITE. 1982. Neurological control of head rotations: electromyographic evidence. J. Neurol. Sci. 55: l-14.