Estimating the time structure of descending activation that generates movements at different speeds

In targeted movements of the hand, descending activation patterns must not only generate muscle activation, but must also adjust spinal reflexes from stabilizing the initial to stabilizing the final postural state. We estimate descending activation patterns that change minimally while generating a targeted movement within a given movement time based on a model of the biomechanics, muscle dynamics, and the stretch reflex. The estimated descending activation patterns predict human movement trajectories quite well. Their temporal structure varies across workspace and with movement speed, from monotonic profiles for slow movements to non-monotonic profiles for fast movements. Descending activation patterns at different speeds thus do not result from a not mere rescaling of invariant templates, but reflect varying needs to compensate for interaction torques and muscle dynamics. The virtual attractor trajectories, on which active muscle torques are zero, lie within reachable workspace and are largely invariant movements when represented in end-effector coordinates. Their temporal structure along movement direction changes from linear ramps to “N-shaped” profiles with increasing movement speed. Author summary What descending patterns of activation drive targeted movements? Based on a model that includes biomechanics, muscles dynamics and the stretch reflex we estimate the descending activation patterns by minimizing their change using optimal control, given a movement target and movement time. The resulting activation patterns predict experimentally observed human movements reasonably well. The temporal structure of the estimated descending activation patterns depends on the speed of the movement, varying from monotonic for slow to non-monotonic for fast movements. This structure reflects the need to compensate for interaction torques and muscle properties. From the model we are able to estimate the virtual attractor trajectories on which all active muscle torques are zero. These lie within reachable workspace, and are relatively uniform across workspace. Their time structure varies from linear ramps for slow movements to N-shaped temporal profiles for fast movements.

Introduction 1 change signal because it also contains a component that reflects the shifted postural 48 state achieved at the end of the movement. This is illustrated in Figure 1 for a shoulder 49 extension movement. Prior to the movement, the flexor muscle is activated only to the 50 extent that it opposes and extensor activation, reflecting co-contraction of flexors and 51 extensors. As the extension movement is initiated, flexor may drop. At the end of the 52 movement, it may return to its initial level of activation (if co-contraction is equivalent 53 in the new posture). A descending "command", which is simply descending activation 54 sent to the motor neurons that drive the flexor muscle, may first drop at the beginning 55 of the extension movement, inducing the muscle's drop in activation. Because the flexor 56 muscle lengthens during the movement, the afferent signals from flexor muscle spindles 57 March 18, 2022 2/23 increase over the movement. Through the reflex loop, this signal would activate the 58 muscle, resisting the movement. Descending activation must drop and remain at a 59 lowered level at the end of the movement, so that it cancels the increased afferent signal 60 and muscle activation may return to its non-movement state. This is solves the 61 movement-posture problem [14]. Behavioral signatures of a persistent component of the 62 descending activation pattern have recently been reported in the force-field adaptation 63 paradigm [15]. 64 The pattern of descending activation is, therefore, not purely change signal, but 65 includes a tonic component that reflects the new postural state. Descending activation 66 itself is not a good candidate for establishing a cost function for movement generation. 67 We propose instead that the cost function depends on the rate of change of descending 68 activation, but definition a change signal. We evaluate that hypothesis by comparing the 69 kinematics predicted by minimizing this cost function with empirical movement data. 70 Our main interest lies not in introducing yet another cost function to describe 71 movement competence, however. Our goal is to use this principle primarily to estimate 72 (or rather: provide a good guess of) the descending activation pattern that is consistent 73 with movement of limbs that are endowed with the stretch reflex. This enables us to 74 investigate how descending activation patterns vary across workspace and across 75 different movement speeds. Given the nonlinear nature of interaction torques, varying 76 movement speed changes their relative contribution during movement. We aim, 77 therefore, to directly observe in the temporal patterns of descending activation 78 signatures of internal models that reflect compensation for interaction torques. 79 This optimal estimate of the temporal patterns of descending activation will also 80 enable us to observe to which extent and over which time interval muscle activation is 81 directly driven by the descending signals. Does the reflex loop contribute to muscle 82 activation during the movement or only once the new postural state has been reached? 83 The stretch reflex defines, in principle, an equilibrium muscle length that determines the 84 postural state at the end of the movement [16]. Does the minimal descending activation 85 pattern define a hypothetical equilibrium configuration of the limb only toward the end 86 of the movement? In the model, we will be able to compute the attractor trajectory [17], 87 on which the limb would move entirely passively as established in experiments in which 88 movement was made passive by a finely attuned external force field [17]. Examining the 89 predicted attractor trajectory across workspace and different movement speeds provides 90 insight into the functional significance of the optimal descending activation patterns.

92
Experimental methods

93
To obtain experimental data for comparison to the theoretical predictions, we performed 94 a standard movement experiment in the manner of [18]. We aimed at naturalistic 95 movements, restraining movement to planar, two-joint action by instruction rather than 96 by mechanically constraining the arm with a manipulandum.

97
Eight movements were chosen to sample a large area of the workspace (Figure 2 The eight examined movements included two movements that covered a distance of 40 117 cm (movement 1 and 2, comparable to movements used in [11]) and six movements that 118 covered a distance of 25 cm. The movements had varying involvement of the two joints. 119 The experiment started with 10 calibration trials, three to estimate joint angles for 120 different arm positions and seven to calibrate the end-effector position for each target. 121 The experimental trials were organized in two blocks, across which two different levels 122 of movement time were imposed (see below). Within each block, movements were 123 randomized. Each block started with a training trial of each of the eight movements, 124 followed by five sessions of 16 trials, two for each movement, for a total of 10 repetitions 125 per movement and movement time. Overall the participants conducted 176 trials.

126
During training trials, participants learned to achieve a target movement time by 127 attending to an auditory metronome. The metronome gave two audio signals of 128 different pitch (750 and 550 Hz). The time interval between the two tones signaled the 129 desired movement time, 400 ms for the fast condition and 800 ms for the slow condition. 130 The two tones were repeated three times, followed by a third tone (620 Hz) that 131 indicated the beginning of the trial. Before the experimental trials, the first two audio 132 signals were played only once before the movement was initiated.

133
Data analysis

134
The motion data was analyzed in Matlab. The trajectories were filtered with a 135 third-order Butterworth filter with a cut-off frequency of 5 Hz. Movement onset was 136 defined as the point in time at which end-effector velocity first exceeded 5% of peak 137 velocity and end-effector acceleration exceeded 5% of maximal acceleration. Movement 138 offset was determined as the point in time when end-effector velocity fell below 5% of 139 peak velocity and end-effector deceleration fell below 5% of peak deceleration.

140
Trials whose movement duration exceeded the median movement duration by more 141 than one third were discarded. If more than four trials for one movement had to be 142 discarded, that movement was deemed unsuccessful and the participant was excluded Theoretical methods 154 We model the reflex control and force generation of a set of six muscles that actuate the 155 shoulder and elbow joints of a planar arm as illustrated in Figure 2. The model largely 156 follows [11], with a few simplifications outlined below. The model takes six descending, 157 time-varying motor commands u desc (t). We determine the "minimal" descending The same biomechanical equations of motion listed above describe the modelled arm 163 dynamics in terms of the joint angle vector, θ = (θ e , θ s ) (elbow, e, and shoulder, s) and 164 its derivatives: The moment of inertia tensor, M, and the Coriolis and Centripetal moments, C, are 166 listed in [18]. Active torques, T , are generated by four monoarticular (musculus 167 pectoralis and deltoid at the shoulder, biceps long head and triceps lateral at the elbow), 168 and two biarticular muscles (biceps short head and the triceps long head). Muscle 169 origins and insertions are taken from [19].

170
Force generation

171
Muscle force generation depends on muscle length, which is determined by the geometry 172 of the arm and the current joint angles. We approximate this dependence linearly, (i = 1, . . . , 6). The parameters, c i (resting length) and c i (moment arm) are taken 174 from [19], and [20], and are listed in Appendix A. We neglect the second order terms of 175 those sources, a reasonable approximation within the studied range of motion (see 176 Figure 1b/c in [21]).

177
Given an activation signal, A i , for each muscle, i, force generation is specified in 178 three steps. An exponential characteristic is weighted by the factor, ρ i , reflecting each muscle's strength that is assumed to vary 180 with the physiological cross-sectional area (PCSA) [19] (values listed in Appendix A).

181
The factor s is constant across all muscles and based on empirical estimates from the 182 March 18, 2022 5/23 cat gastrocnemius muscle [22]. In a second step, calcium kinetics is modeled as a 183 critically damped second order low pass filter [11] 184 τ 2M + 2τṀ + M =M (6) with a single parameter, τ , listed in Table 5. The third step takes into consideration the 185 dependence of muscle force, F i , on the muscle contraction rate [23] and passive muscle 186 elasticity, where c i is the passive resting length of the muscle and the factors f 1 to f 4 and k were 188 adjusted to reproduce the results of [11] (values listed in Appendix A).

189
The active joint torques, T , are obtained from the muscle forces, F i , taking the 190 moment arms into account: The activation, A i , of each muscle is assumed to reflect the Descending activation, 193 u desc,i acts as the threshold of a reflex loop that controls muscle activation, A i , for each 194 muscle, i [11]: where [·] + is a semi-linear threshold function. The parameters γ and µ reflect the 196 sensitivity of muscle spindles to muscle length and its rate of change. The resting level 197 of the α-motorneuron in the absence of descending activation is h i < 0. Note that at 198 rest (l i = 0), descending activation (more precisely, −(u desc,i + h)/γ) can be construed 199 as a threshold length of the muscle as is done in the λ-model of posture [16]. When under the boundary conditions To guarantee the smoothness of the modeled movement, 216 we also require that initial and terminal velocity and acceleration are zero:  Table 1. Movement time (T ) and its standard deviation (T SD ) across participants for each of the eight movements in the two movement conditions, slow vs. fast.
upper and lower bounds for joint angles and descending commands, and thus for muscle 219 forces, in order to stay within biomechanical and physiological limits.

220
The optimization problem was solved numerically using the control vector 221 parameterization approach [24,25] as described in the Appendix C. In the model, it is possible to directly compute the virtual attractor trajectories that are 236 accessible empirically only through sophisticated instrumentation [17]. At any moment 237 in time, the virtual attractor state is the set of joint angles and associated muscle 238 lengths, at which the active torques contributed by all muscles sum to zero. In the 239 model, we compute this state at each point in time given a current descending 240 activation vector, u(t), and current rates of change of all muscle lengths, d l(t)/dt.

241
Virtual attractor trajectories were obtained in joint space, but then transformed into 242 hand position space using the kinematic model.

256
Hand paths from experiment match those from the model qualitatively ( Figure 3).

257
For all movement conditions, the small deviation from a straight path was consistent in 258 direction and amount of curvature between model and experiment. Because participants 259 have different segment lengths and were not, therefore, identically positioned relative to 260 the targets, assembling mean paths required alignment across participants. This was 261 done by centering all hand paths in the shoulder position. As a result, experimental 262 paths are shifted relative to the simulated paths.

263
The two-dimensional hand trajectories through space are shown in Figure 4 for the 264 the two matching movement times conditions of experiment and model. Normalizing  Generally, the model of the stretch reflex, Equation 9, implies that when a muscle 347 lengthens, the same amount of descending activation generates more muscle activation 348 than when a muscle shortens. Because acceleration is caused by muscles that shorten 349 and breaking is caused by muscles that lengthen, this implies that more change of 350 descending activation is needed to accelerate than to break a movement. Minimizing 351 descending activation change thus leads to slower/lower acceleration than deceleration. 352 Muscle physiology would tend to compensate for this effect as contracting muscles 353 produce more force than lengthening muscles, reducing the amount of descending 354 activation needed to bring about acceleration. Hill type muscle models capture this 355 physiological effect [20,28], and would thus reduce or potentially reverse this discrepancy. 356 In contrast, the simplified model of muscle force generation that we used [11], amplifies 357 the discrepancy because in this model, the same level of muscle activation generates 358 more force for lengthening than for shortening muscles. We chose this simplified model 359 to enable comparison to that earlier work (see below). The simplified model also 360 sidesteps the additional redundancy implied in Hill-type models by the non-unique 361 mapping between the length of the contractile element and the joint angle.

362
What does the cost function mean?

363
The capacity of the model to predict the observed kinematics and its invariance across 364 rate and workspace is remarkable given that no physical or physiological factor enters 365 into the cost function. This is in contrast to other approaches, such as minimal torque 366 change [29] and minimum effort [13]. The latter model in effect minimizes muscle 367 activation. The asymmetry of the velocity profiles is not consistently captured by that 368 model either (see their Figure 4b, which is not directly compared to experiment, 369 however).

370
Thus, merely minimizing the temporal complexity of the descending activation 371 pattern captures some of the regularity of human movement kinematics. This is 372 consistent with the general theme that the peripheral reflex circuitry and muscle 373 dynamics simplify motor control [10,30].

374
This theme was at the center of the modeling study from which we adopted the 375 reflex and simplified muscle model [11]. That study showed that simple "ramp" 376 descending activation patterns (framed by those others as virtual trajectories) are 377 sufficient to generate realistic movement patterns and stiffness profiles. Our approach 378 directly speaks to that observation. In fact, the optimal control problem as framed here 379 leads to such ramp-like solutions in the limit case of a trivial plant. Interpreting the 380 descending activation as threshold muscle length, λ [16], the trivial plant is realized 381 when muscle lengths simply track threshold lengths, l(t) = λ(t). In that case, the 382 optimal control problem can be easily solved analytically and leads to a linear ramp in 383 λ(t) that connects the initial to the final value for each muscle. Since we initialized 384 optimizations with exactly this linear ramp, the resulting descending activation patterns 385 results may indeed be interpreted as the solutions closest to the linear ramps while 386 compatible with the imposed movement target and time. These solutions deviate from 387 the linear ramp because the real plant deviates from the trivial one due to biomechanics, 388 muscle dynamics, and reflex delays.

389
The ramp-like profiles we found for slow movements approximates the limit case of 390 the trivial plant, although the ramp ends before the movement ends. This is consistent 391 with experimental estimates of the time when the descending signal reaches its terminal 392 level [31]. Those experiments actually involved fast movements, and did not estimate the 393 temporal profile of the descending activation pattern. Our estimates for fast movements 394 identified temporally more complex, non-monotonic descending activation patterns.

395
Earlier work hypothesizing ramp-like profiles for the same neuro-muscular model used 396 10/23 here [11] did not explore movement speeds and portions of work space at which we 397 found these non-monotonic estimates. Thus, the modeling work reported here clarifies 398 earlier debates about the complexity of descending signals required to bring about 399 movement.We do not, however, seriously propose that the nervous system minimizes the 400 temporal complexity of descending activation patterns. The optimal control approach 401 was used primarily as a tool to provide reasonable estimates of the descending 402 activation patterns, which we may then analyze for their regularities and properties.

403
The temporal structure of descending activation patterns 404 We discovered a range of different time structures of the minimally changing descending 405 activation patterns (Fig. 7, 8). deceleration of the movement acquire this more strongly modulated temporal structure 413 for faster movements, consistent with the need to overcome inertial and interaction 414 torques ( Fig. 7 and 8). Qualitatively, these descending activation patterns resemble the 415 "N-shapes" postulated in [27] (although that earlier work entailed single joints and an 416 inadequate linear model as argued in [32]). Descending activation patterns of muscle 417 pairs that contribute weakly to the acceleration and deceleration of motion tend to be 418 strongly modulated in time.

419
The contribution of biarticular muscles that span the shoulder and elbow joints  co-contraction emerged from the optimization. By varying the imposed initial/final level 436 of co-contraction we found that these levels had little influence on the descending 437 activation patterns generated by the optimization and on the predicted movement.

438
The observed difference between temporal profiles of descending activation patterns 439 of slow versus fast movements refutes the idea that movement generation can be 440 explained by a single movement template or primitive, that would be scaled to the 441 desired movement time. Figures 7 and 8) already use normalized time so a simple 442 scaling would imply identical temporal structure for different movement speeds. When 443 descending activation patterns are specifically shaped to anticipate velocity dependent 444 torques and muscle properties, then the neural processes generating those temporal 445 profiles could be considered to reflect an "internal model" as postulated in 446 computational theories of motor control [6,33].

447
The virtual attractor trajectory 448 The time course of the joint angles and associated muscle lengths which all active 449 torques generated by the muscles are zero is called the virtual attractor trajectory 450 by [17]. Unlike in experiment, it is easy to determine these virtual attractor trajectories 451 in the model. We found these to be astonishingly uniform across the different 452 movements when plotted in end-effector space in coordinate systems that are aligned 453 with each movement (Figure 9). Along the movement, the virtual attractor trajectories 454 of slow movements are linear ramps that reach all the way to the end-point of the

461
It is important to note that all virtual attractor trajectories lie within the reachable 462 work space. Thus, at every moment in time during a movement, a realisable muscle 463 length exists for all muscles at which descending activation is compensated for by 464 activation from the stretch reflex (neglecting co-contraction). If the hand was to exactly 465 trace the virtual trajectory (as brought about in the experiments of [17] by applying 466 exactly the right external torques), activation that descends would be exactly equal to 467 the activation that the reflex loop contributes (again, except for co-contraction). In a 468 sense, the distance between the real and the virtual hand trajectory reflects the 469 difference between these two contributions. The observed patterns suggest that the 470 contribution of the reflex to movement generation is of the same order of magnitude as 471 the contribution of descending activation itself, especially for slow movements.

472
It may be tempting to think of the invariant shape of the virtual trajectory as a 473 simple control strategy. Note, however, that the descending activation patterns that the 474 brain's neural networks must generate are less invariant than the virtual trajectories 475 that reflect the contributions of biomechanics, muscle dynamics, and reflex loops. The 476 invariance of the virtual trajectory thus essentially is the dynamic reflection of the 477 invariance of human movement kinematics across work space and speed.

478
Limitations 479 Clearly, this model provides only a first rough approach to using models of the periphery 480 to estimate descending activation patterns. Using more realistic and detailed muscle 481 models and including more complete accounts of spinal reflexes is a future challenge [10]. 482 Muscle redundancy was also only addressed summarily, but minimizing the change 483 of activation to redundant muscles. A more thorough-going analysis of muscle 484 redundancy may require new experimental approaches to enable estimation of the 485 descending activation of individual muscles. Similarly, kinematic redundancy [34] adds 486 an important next level of complexity that would need to be addressed by the approach 487 we sketched here.

488
For computational reasons, we neglected delays in the reflex loop. To assess the 489 error induced, we input the estimated descending activation patterns into a model that 490 included typical sensory delays (as quantified in [11]). We found that the resultant 491 kinematics differ little from the kinematics of the model without delay. This is likely 492 due to the fact that in the most dynamic portions of the movement, early acceleration 493 and late deceleration, joints and muscle lengths vary slowly.   Table 2. Moment arms and resting lengths of the muscles are taken from [20]. Muscle lengths are taken from [19]. 3.87 Table 4. Values for the physiological cross section areas for the six modeled muscles are taken from [19].  Table 5. Values of parameters that specific to the neuro-muscular model were taken from [11] and adjusted such that the results of [11] could be reproduced. Parameters f 1 , f 2 , f 3 f 4 and c are fitting parameters for the active and viscous part of force generation. The passive stiffness is scaled with k. τ is the time constant of the differential equation for the muscle force generation.

B Co-contraction 529
The amount of co-contraction, C, specifies how much all threshold lengths, λ, must be 530 shifted such that the forces generated in equilibrium by individual muscles reach that 531 given level without changing net joint torques. We separate the descending signals λ 532 into a reciprocal part λ r and a co-contraction part λ c with λ r + λ c = λ. Furthermore, 533 we consider the stationary state without external torques. Under those conditions, the 534 reciprocal part of the reference command equals the actual muscle lengths, λ r = l, and 535 does not contribute to the generation of forces. To determine the co-contraction 536 component of the reference command, λ c , we compute the threshold lengths at which 537 the force generated by each muscle is increased by an average of C (in Newton) by 538 solving: Here, T i,e and T i,s are the elbow and shoulder torques induced by the i-th muscle.

540
C Control vector parameterization 541 We used the control vector parameterization approach [24] to optimize the time course 542 of the reference commands λ(t). Each of the six continuous reference commands, λ i , 543 was sampled at N = 20 moments in time, leading to a total number of 120 parameters 544 to optimize. Optimization was performed with respect to the objective function defined 545 in section using the iterative nonlinear programming solver fmincon. Initialization was 546 chosen so that the reference commands sampled a linear ramp from the initial to the 547 target position in the Cartesian plane. After each optimization step, cubic splines were 548 used to interpolate the reference commands between the discrete sample points to    . In a simple model of the stretch reflex (illustrated on the bottom right), activation that descends from the brain is combined with length and velocity feedback sensory signals from muscle spindles. When these summed inputs exceed threshold, muscle activation is induced. During shoulder extension, flexor muscle activation first drops from an initial level that reflects the level of co-contraction of flexor and extensor muscles, then rises, contributing to the breaking of the movement, and finally returns to its initial value. Because the flexor lengthens, the feedback signal is larger at the end of the movement than initially. The descending activation pattern is lowered by a matching amount, so that feedback and descending activation cancel, reducing muscle activation and suppressing the flexor's resistance to the new posture. All curves are from the model, which finds the descending activation pattern that changes minimally to achieve the new posture within a given movement time (see Results Section). Top: Eight movement conditions (the labels 1-8 are positioned near the end-point of each movement) sample the space of planar movements performed with two degrees of freedom, shoulder and elbow. Bottom: The model includes shoulder and elbow joints in the horizontal plane, articulated by two mono-articular elbow joint muscles (red), two mono-articular shoulder muscles (blue) and two bi-articular muscles that span elbow and shoulder joint (green). The illustrated geometry leaves muscle lever arms invariant across work space.    Simulated and experimental active joint torques across normalized time for all eight movements (rows) for the shoulder (green) and elbow (lilac). The first pair columns show the slow condition, the second pair of columns shows the fast condition. In each pair, the simulation comes first, the experiment second. The shaded areas in the plots of the experimental data shows the standard deviation.