The cost of being stable: Trade-offs between effort and stability across a landscape of redundant motor solutions

Current musculoskeletal modeling approaches cannot account for variability in muscle activation patterns seen across individuals, who may differ in motor experience, motor training, or neurological health. While musculoskeletal simulations typically select muscle activation patterns that minimize muscular effort, and generate unstable limb dynamics, a few studies have shown that maximum-effort solutions can improve limb stability. Although humans and animals likely adopt solutions between these two extremes, we lack principled methods to explore how effort and stability shape how muscle activation patterns differ across individuals. Here we characterized trade-offs between muscular effort and limb stability in selecting muscle activation patterns for an isometric force generation task in a musculoskeletal model of the cat hindlimb. We define effort as the sum of squared activation across all muscles, and limb stability by the maximum real part of the eigenvalues of the linearized musculoskeletal system dynamics, with more negative values being more stable. Surprisingly, stability increased rapidly with only small increases in effort from the minimum-effort solution, suggesting that very small amounts of muscle coactivation are beneficial for postural stability. Further, effort beyond 40% of the maximum possible effort did not confer further increases in stability. We also found multiple muscle activation patterns with equivalent effort and stability, which could underlie variability observed across individuals with similar motor ability. Trade-off between muscle effort and limb stability could underlie diversity in muscle activation patterns observed across individuals, disease, learning, and rehabilitation. Author summary Current computational musculoskeletal models select muscle activation patterns that minimize the amount of muscle activity used to generate a movement, creating unstable limb dynamics. However, experimentally, muscle activation patterns with various level of co-activation are observed for performing the same task both within and across individuals that likely help to stabilize the limb. Here we show that a trade-off between muscular effort and limb stability across the wide range of possible muscle activation patterns for a motor task could explain the diversity of muscle activation patterns seen across individuals, disease, learning and rehabilitation. Increased muscle activity is necessary to stabilize the limb, but could also limit the ability to learn new muscle activation pattern, potentially providing a mechanism to explain individual-specific muscle coordination patterns in health and disease. Finally, we provide a straightforward method for improving the physiological relevance of muscle activation pattern and musculoskeletal stability in simulations.


Author summary
Current computational musculoskeletal models select muscle activation patterns that minimize 24 the amount of muscle activity used to generate a movement, creating unstable limb dynamics.

25
However, experimentally, muscle activation patterns with various level of co-activation are 26 observed for performing the same task both within and across individuals that likely help to 27 stabilize the limb. Here we show that a trade-off between muscular effort and limb stability across 28 the wide range of possible muscle activation patterns for a motor task could explain the diversity 29 of muscle activation patterns seen across individuals, disease, learning and rehabilitation.

30
Increased muscle activity is necessary to stabilize the limb, but could also limit the ability to learn 31 new muscle activation pattern, potentially providing a mechanism to explain individual-specific 32 muscle coordination patterns in health and disease. Finally, we provide a straightforward method

Introduction 36
Little is known about the landscape of redundant motor solutions for sensorimotor tasks, hindering 37 our ability to interpret variability in muscle activation patterns seen across individuals, who may 38 differ in experience, training, or neurological health. In general, we refer to motor solutions as the 39 spatiotemporal pattern of motor neuron activity that are reflected in electromyographic (EMG) 40 signals during posture and movement. More specifically, we have demonstrated that individual-41 specific motor modules-defined as library of fixed spatial patterns of muscle activity that execute 42 a motor task-can vary considerably across individuals. Experimental EMG recordings reveal 43 variations across individuals in the degree to which muscles are co-activated in motor modules 44 that perform similar tasks [1,2]. Further, differences in motor module structure are amplified when 45 examining individual across the spectrum of motor ability: long-term motor training, such as in 46 ballet dancers, tends to decrease the number of muscles used to execute a given task [3], whereas neurological impairments tend to increase the co-activation of muscles used to perform 48 gait and balance tasks [4][5][6][7][8]. For reactive balance tasks in cats, we explicitly showed individual-49 specific differences in motor modules that receive similar neural commands and produce 50 equivalent motor outputs in healthy, trained, animals [9]. Through analysis of musculoskeletal 51 models, we found the biomechanics of the limb and task requirements to impose few constraints 52 on feasible muscle activation patterns [10][11][12][13]. These results indicate that there is ample room in 53 the nullspace, i.e. set of possible motor solutions, from which individuals may employ motor 54 modules that deviate from the muscle coordination pattern predicted by single optimality criteria 55 such as minimizing muscle effort [14][15][16]. However, we understand little about what might drive 56 the selection of particular individual-specific patterns within the nullspace of a given task.

58
While minimizing muscular effort has been established as a principle driving muscle 59 activation patterns for movement, increasing muscle effort to stabilize the limb is also an important 86 a musculoskeletal model of that cat hindlimb, we demonstrated that limb stability can be 87 modulated by the selection of diverse patterns of muscle activation pattern within the nullspace 88 of equivalent solutions in terms of force output [40]. Stability was evaluated by computing the 89 eigenvalues of the linearized system dynamics, incorporating both rigid body dynamics and force-90 length and force-velocity properties of activated muscles [36,40]. Our prior work only examined 91 stability of randomly sampled solutions throughout the nullspace, without consideration for 92 muscular effort. To our knowledge, no study has characterized the functional landscape of stability 93 along the trajectory between minimum-and maximum-effort solutions where effort level is 94 systematically varied. Here our goal was to quantitatively express trade-offs between effort and 95 stability to understand how motor solution might differ across tasks and individuals.

97
Given that individuals may seek to modulate muscular effort and/or limb stability [41,42],

98
it is also not known whether one could learn a better motor solution by making small changes to 99 a muscle activation pattern within the nullspace of motor solutions. The observation that 100 individuals maintain characteristic motor solutions, e.g. habitual [43][44][45] or pathological [46][47][48],

101
suggests that it may be difficult to find nearby solutions that improve the functional properties 102 while also satisfying task requirements. In particular, without a priori knowledge of what direction 103 of change in muscle space would satisfy the task constraint, how easy is it to find a neighboring 104 solution that satisfies the tasks constraint and desired changes in the functional properties?

105
Describing the local landscape around particular motor solutions with respect to effort and stability

129
In summary, our overall objective was to explore the nullspace of motor solutions for a 130 quasi-static force production task ( Fig 1A) with respect to the two functional properties of effort 131 and stability using a 3D musculoskeletal model of the cat hindlimb ( Fig 1B) 2 at the knee, and 2 at the ankle) and 31 Hill-type muscle actuators (see Table 1

159
Trade-off between effort and stability along a continuum of solutions in the 160 nullspace 161 We found dramatic differences in the effort and stability of the minimum-and maximum-effort 162 solutions for generating an extensor force typically observed during quiet standing as well as 163 during postural responses to perturbation in cats ( Fig 1A). When normalized as percent of the 164 maximum possible effort level (i.e. 100%), the same task could be achieved with the minimum-165 effort solution at only 8.80% effort of the maximum-effort solution. We then evaluated the stability 166 of each solution using a stability metric defined as the maximum real part of the eigenvalue of the We found a reciprocal relationship between effort and stability when stability was 194 examined along the continuum of solutions between the two the extremes identified by the 195 minimum-and maximum-effort solutions (Fig 2A). We identified the set of muscle activation 196 patterns along a null-path, defined as the linear combinations of minimum-effort and maximum-197 effort solutions ( Fig 2B) for producing the same endpoint force vector. As such, the minimum-

208
The slope of the null-path indicates the sensitivity of the system's postural stability to 209 changes in muscular effort, which is very steep for low effort solutions and insensitive for solution 210 of greater than ~40% effort (Fig 2A). The most sensitive region was between the minimum-effort 211 solution and the first stable solution, where small increase in effort resulted in rapid improvement 212 in stability (e.g. Min-E in Fig 2C). From the most stable solution to the maximum-effort solution, 213 however, stability was almost insensitive to change in effort (e.g. Most-stable in Fig 2C). More 214 importantly, beyond the most stable solution (E=38.6%), further increases in effort decreased 215 stability (e.g. Max-E in Fig 2C).

217
Local landscape of effort and stability vary along the null-path 218 Because the null-path only identifies solutions along one particular direction in the high 219 dimensional nullspace (24, for our musculoskeletal model), we further explored the local 220 landscape around 51 solutions along the null-path (Fig 3). We were particularly interested in  normally distributed (red) about the activation of the seed solution ('x) with step size (=1σ) 239 corresponding to 5, 10, and 25% of the feasible range were generated. Neighboring solutions 240 (dots) were then found by projecting theses perturbed patterns to the solution manifold. Resulting 241 distributions are shown in thick solid line for 5% step size, thin solid line for 10% step size, and 242 dotted line for 25% steps size. (C) Two-dimensional 95% confidence ellipses were found for each 243 set of neighboring solutions, i.e., for each seed solution (four representative examples in 'x') and 244 steps size (5%, 10%, and 25% in thick solid, thin solid, and dotted line, respectively). The ellipse 245 was centered at the mean effort and stability, with major axes in the direction of maximum 246 covariance of neighboring solutions in the functional property space. The local landscape 247 explored by the neighboring solutions closely followed the null-path (gray line). Changes in 248 functional property occurred mainly in stability near low effort region, e.g. the minimum-effort ('x' 249 in light blue) solution, while it was very difficult to change stability near high effort region, e.g. the 250 maximum-effort solution ('x' in black). Actual neighboring solutions are shown on the right panels 251 (note different scales), with distributions shown in thick solid, thin solid, and dotted lines for 5%, 252 10%, and 25% step size, respectively. [a.u.] stands for arbitrary unit. (D) With increased step size, 253 the orientation of the 95% ellipse deviated more from the slope of the original seed solution on 254 the null path (ΔSlope), and larger region (Area) with more different effort and stability (Offset) 255 were found. Explored solutions also spanned wider range of directions around the seed solution 256 (Eccentricity) with increased step size. 257

258
Local searches around seed solutions on the null-path preserved the functional properties

277
As step size increased, solutions differed more compared to the seed solution ( Fig 3D).

278
Across local searches from all 51 seed solutions, the difference between the major axis of each 279 ellipse and the slope of the null-path for each seed solution increased with step size (ΔSlope; 280 χ 2 (2)=12.4, p=0.002). In particular, the change in orientation of the ellipses from the slope of the 281 null-path were greater with 25% step size, compared to 5% (p=0.002) and 10% (p=0.041). As 282 step size increased, the area (χ 2 (2)=76.1, p=3.05e-17) of the ellipse and the offset of the center 283 from the seed (χ 2 (2)=76.8, p=2.09e-17) also increased, whereas eccentricity decreased 284 (χ 2 (2)=6.85, p=0.033), indicating that solutions that are further from the null-path in the functional

285
property space were explored. The area for the 25% step size search was greater than for the 5%

292
Multiple local solutions can have equivalent functional properties 293 We found multiple solutions that had very similar effort and limb stability near the most- hip extensors (e.g. gmed and gmin, see Table 1 for abbreviations) differed between the two 301 solutions (Fig 4B, top). Hamstring muscles (e.g. sm, st) were highly activated in one solution (Fig  4B, magenta), which resulted in activation of hip flexor/knee extensor (e.g. rf) for torque balance 303 at the knee. However, in the other solution ( Fig 4B, red), vm was mainly used for knee extension.

304
Large activation of mg and lg, along with co-activation of sol and ta, were present in one solution 305 (Fig 4B, magenta), whereas the other solution ( Fig 4B, red)   with a heuristic search, and using all solutions explored in this study (area encompassed in 358 orange). The shape of the Pareto front closely followed the null-path (gray), but was consisted of 359 solutions that were more stable than the null-path. (B) Muscle parameters were varied to 360 investigate how other intrinsic stabilizing mechanism and altered condition in force-length 361 relationship affect the mapping in the functional property space. Compared to the generic model 362 with muscles operating at 65% normalized fiber length (solid gray), adding short-range stiffness 363 (dotted gray) drastically improved stability for all solutions on the null-path. Solutions on the null-364 path for the model with muscles operating at 95% normalized fiber length i.e., plateau region of 365 force-length curve, was always unstable (purple), which could not be stabilized even with short-366 range stiffness (dotted purple). 367 368 Effects of muscle parameter variations on functional properties 369 Altering intrinsic muscle properties changes the limb stability but not the qualitative effort-stability 370 trade-off relationship (Fig 5B). Stability drastically improved when intrinsic stabilizing mechanism 371 such as short-range stiffness was added to the generic model with fibers operating at 65% 372 normalized length (Fig 5B, dotted  observed during gait and balance tasks in healthy humans and animals [1,2,9,60,61,105].

515
Rather than being at an absolute optimal solution in terms of effort, the additional criteria of local 516 stability suggest that individuals use "slop-timal" [106,107], or good-enough [53] motor solutions.

519
Musculoskeletal model and target endpoint force 520 We used a detailed musculoskeletal model of a cat hindlimb [108] to examine muscle activation 521 patterns to produce an experimentally observed endpoint force during postural response in a cat

549
For the target endpoint force vector, we used an extensor force vector ( Fig 1A) that was

Metrics for effort and stability 556
In order to evaluate functional properties of redundant muscle patterns that generate the same 557 endpoint force and to map a functional property space in terms of effort and stability, we defined 558 quantitative metrics for each criterion. We defined the metric for effort (E) to be the square-root of 559 sum of squared activations (Eq 3), which is equivalent to summing muscle stress [15, 68]:

560
(3). mapping (Eq 2), using quadratic programming. The effort for any muscle activation pattern examined in this study was then normalized to percent of the global maximum, i.e., E of the 566 maximum-effort solution.

568
We defined the metric for stability using Lyapunov stability of the linearized model. The 569 full nonlinear system (Eq 1) was linearized about a static equilibrium point, defined by a muscle 570 activation pattern that satisfies the endpoint force generation (Eq 2), using software 571 Neuromechanic [114]. Specifically, the system equation incorporated joint torques generated by 572 muscles: The system was linearized by numerically computing the partial derivatives with respect to 575 kinematic states using first-order Taylor-series expansion to obtain the state space representation: The state matrix (A) was used to calculate the eigenvalues ( ) of the linearized system. For a  578 given muscle activation pattern, the metric for stability (S) was defined as the maximum real part corresponding to the 6 eigenvalues that are near zero (similar to rigid-body modes) which are very small in magnitude. These 6 eigenvalues typically were less than 1e-5 in magnitude, and 590 thus were not relevant in physiological time scales: the time at which the magnitude of a perturbed 591 response is reduced to 50% is longer than 6.9e4 seconds. In contrast, the 8 eigenvalues that 592 were considered in S for a given solution were typically larger than 1e-3 in magnitude.

594
A solution of the system defined by a given muscle activation pattern was determined 595 "stable" if S<0, and "unstable" if S>0. Further, because the magnitude of S predicts the rate at 596 which a perturbed system will return to (if S<0), or deviate from (if S>0) the equilibrium, a solution 597 is to be "less stable" for greater value of S, and "more stable" for smaller value of S. This metric 598 from system theory, i.e., Lyapunov indirect or linearization method, has been shown to predict the 599 behavior of perturbed nonlinear systems in simulations [36,40].

612
we defined the unique path between the two solutions in the muscle activation space as the null-613 path.

615
We mapped the null-path onto the functional property space by evaluating the effort (E) 616 and stability (S) of each intermediate solution, revealing the projection from the linear muscle 617 activation space to the functional property space. In particular, we were interested in whether the 618 minimum-effort solution, representing the least amount of co-activation for a given task, is 619 unstable, and whether the maximum-effort solution, with highest level of co-activation, is stable.

620
Further, this null-path was used to demonstrate whether an unstable solution can be made stable 621 by following a defined direction in muscle activation space, and to quantify the amount of effort 622 necessary to make the minimum-effort solution stable.

624
To quantify the sensitivity of stability to change in effort, we calculated instantaneous slope The step size was defined for each muscle as a percentage of the feasible muscle 641 activation range (Sohn et al., 2013) for the target endpoint force production, so that the amount 642 of change in any muscle is normalized. To vary the extent to which area in the functional property 643 space is explored by the neighboring solutions around a given seed solution, we used step sizes 644 of 5, 10, and 25% of the feasible muscle activation range.

646
We generated 262 neighboring solutions ( ) for each seed and step size using the neighbor e  647 following four steps. In summary, our goal was to explore the functional property space in every 648 possible direction around a given seed, reaching the maximum extent in distance for a given 649 amount of changes allowed in the muscle activation space. To this end, we generated sets of 650 solutions that were randomly distributed around a given seed with varying distance in muscle

669
Note that we used a scaling matrix Q, which is a diagonal matrix that penalizes the difference 670 between perturbed and projected muscle activation. Elements of Q were weighted inversely 671 to the feasible range to prevent projection only occurring in muscles with small feasible ranges.

672
Each optimization was subject to an equality constraint for satisfying the torque requirement: for projected solution (Eq 11-13), where the problem is reduced to solving for a 30-dimensional 695 muscle activation vector that minimizes the distance to either lower or upper limit vector (Eq 696 13), specified for a given step size. Although solutions generated this way were not 697 guaranteed to satisfy the torque requirement, especially when small step sizes were allowed, 698 more than 57 (60±1.1 across seed and step size) limit solutions were generated for each seed 699 and added to the set of neighboring solutions of the seed for a given step size.

701
Characterization of local landscape 702 We characterized the local landscape around each seed solution using the neighboring solutions 703 mapped in the functional property space. We were particularly interested in whether 704 characteristics of the original seed solution on the null-path are preserved as step size increased.

705
To this end, we identified a 95% confidence ellipse in the two dimensional space of effort and 706 stability (Fig 3C), using the covariance matrix of effort (E) and stability (S) values evaluated for 707 each set of neighboring solutions, i.e., for each seed and each step size. The eigenvector of the covariance matrix was used to define the major axis of the ellipse, corresponding to the direction 709 in which effort and stability varied maximally, and the minor axis that is orthogonal to the major 710 axis. Note that such procedure is equivalent to performing a principal component analysis. The

711
center of the ellipse was defined at the mean effort and stability values for given set of neighboring 712 solutions. The radii of major and minor axes were defined by scaling the maximum and minimum 713 eigenvalues of the covariance matrix, respectively, with chis-squared value for 95% percentile.

714
Note that while such ellipse may not encompass or be fully occupied by the actual neighboring properties, we compared multiple solutions that are within a small range of effort and stability (Fig   745  4). Specifically, we were interested in whether near-maximal stability could be achieved using 746 vastly different muscle patterns. Thus, among all neighboring solutions investigated, we examined 747 solutions that were nearest to the most-stable solution (Fig 4A, 'x' in dark blue) on the null-path 748 between the minimum-and maximum-effort solutions. In order to avoid selecting solutions that 749 are inherently similar in muscle activation pattern, we first sorted 15 solutions that were nearest 750 in Euclidean distance to the local maximum stability solution from each set of neighboring 751 solutions around a given seed, across step sizes. Among these 765 solutions (15 from each of 752 51 seeds), we selected solutions that had effort and stability within 10% and 1% difference, 753 respectively, from the most-stable solution on the null-path. This resulted in total 166 solutions 754 near the most stable solution on the null path ( Fig 4A, gray dots).
755 To quantify the spatial amount of redundancy, we computed percentage of range spanned 757 by the above nearest solutions for each muscle with respect to corresponding feasible muscle 758 activation range. To examine the redundancy in the muscle activation space, we also computed 759 the dimensionality within the activation set [40]. Dimensionality was computed using principal 760 component analysis, defined as the number of principle components required to explain over 95% 761 variance in data (i.e., a 31x166 matrix).

763
To determine whether similar level of stability was achieved by different types of stabilizing principle components required to explain over 95% variance in data (i.e., a 14x166 matrix).

774
Global maximum-stability and Pareto front 775 We found a global maximum-stability and defined the edge that connected the minimum-effort 776 solution and global maximum-stability solution as the Pareto front and compared the values and 777 shape to the null-path between the minimum-and maximum-effort solutions. To examine an 778 explicit trade-off for minimizing effort while maximizing stability, we identified Pareto front that 779 quantifies the maximum level of stability that can be achieved for a given amount of effort, or vice 780 versa, the minimum amount of effort required to achieve certain level of stability. To define this Effects of muscle parameter variations on functional properties 792 We further investigated how altered condition in force-length relationship and other possible 793 intrinsic mechanisms that contribute to stability affect the mapping in the functional property space.

794
First, stability of the original 51 solutions on the null-path were evaluated in a model implemented 795 with a short-range stiffness in muscles [83,117] to represent the instantaneous force-length 796 behavior. However, the short-range stiffness model was modified because we assumed tendons 797 to be inelastic in the static mapping (Eq 2). Hence, our short-range stiffness model represented 798 the instantaneous stiffness only from the active muscle, assuming an infinite stiffness in the 799 tendons that were in series. In result, we used a force-relationship curve with a constant slope 800 [83]. In addition, we did not neglect the force-velocity relationship that provided damping as in the 801 generic model because its contribution may be large in the linearization process where we used 802 numerical perturbation method.

804
Secondly, we computed another set of 51 solutions that constituted a null-path but with a 805 modified model where fiber lengths of all muscles were set at 95% optimal fiber length, i.e.,

806
plateau region of the force-length relationship curve [112] and thus more prone to destabilization.

807
Note that these solutions are not the same as the original solutions, because the scaling matrix 808 for isometric force generation ( ) in Eq 2 is different. However, each of the 51 solutions in the AFL F 809 two models (65% and 95%) were very similar (cosine of the angle between the two muscle 810 activation vectors was 0.998±0.002 across 51 solution pairs). In addition, we also evaluated the 811 null path of the 95% model with added short-range stiffness properties of active muscles using 812 the same method described above.