Optimal Workloop Energetics of Muscle-Actuated Systems: An Impedance Matching View

Integrative approaches to studying the coupled dynamics of skeletal muscles with their loads while under neural control have focused largely on questions pertaining to the postural and dynamical stability of animals and humans. Prior studies have focused on how the central nervous system actively modulates muscle mechanical impedance to generate and stabilize motion and posture. However, the question of whether muscle impedance properties can be neurally modulated to create favorable mechanical energetics, particularly in the context of periodic tasks, remains open. Through muscle stiffness tuning, we hypothesize that a pair of antagonist muscles acting against a common load may produce significantly more power synergistically than individually when impedance matching conditions are met between muscle and load. Since neurally modulated muscle stiffness contributes to the coupled muscle-load stiffness, we further anticipate that power-optimal oscillation frequencies will occur at frequencies greater than the natural frequency of the load. These hypotheses were evaluated computationally by applying optimal control methods to a bilinear muscle model, and also evaluated through in vitro measurements on frog Plantaris longus muscles acting individually and in pairs upon a mass-spring-damper load. We find a 7-fold increase in mechanical power when antagonist muscles act synergistically compared to individually at a frequency higher than the load natural frequency. These observed behaviors are interpreted in the context of resonance tuning and the engineering notion of impedance matching. These findings suggest that the central nervous system can adopt strategies to harness inherent muscle impedance in relation to external loads to attain favorable mechanical energetics.


Supporting Material Bilinear Model for Muscle Force
In the Materials and Methods section, we assumed that for the oscillatory motions considered, muscle mechanical impedance can be approximated by a bilinear model. The bilinear model described by Equation (3), is repeated here: How did we arrive at this approximation, and what is the accuracy of this model? Here we provide a brief summary based on [1].
To arrive at a muscle impedance model, we conducted measurements to characterize the mechanical impedance of the Plantarus longus muscles, similar to the ones used in our investigation (but not the particular muscles used to test the hypothesis because of the limited fatigue life of the muscles). We used the same materials and methods described in the main text. This includes muscle type, harvesting procedure, and the use of the servo-mechanical muscle testing apparatus. There were two differences, however, due to the fact that measurements focused on uncovering the dependence of muscle force on mechanical state, rather than estimating workloop energetics. First, sinusoidal motions were imposed on the muscles, and electrical stimulation was delivered at predetermined points in time. Thus, the experimental setup was akin to Figure 1A in the main text. The electrical stimulus points were not synchronized with the oscillations, but rather were chosen to span various displacement-velocity combinations along the sinusoidal trajectory of the muscle. Thus, the timing of electrical stimulation trigger points was uncorrelated with the imposed motion. Sample oscillatory motions are illustrated in Figure  S1A. The second difference is that the electrical stimulus parameters were fixed, and were not dependent on prior optimal control computations. The electrical stimulus parameters were set to at duration of 50 msec, consisting of 10 spikes at 200 Hz, with a spike duration of 100 micro-seconds.
To explore the relative importance of various impedance factors, we used the acquired data to fit a general impedance function in the form of: While additional higher-order model terms could have been added, this general model can capture many facets of the impedance characteristics of the muscle, including passive stiffness ( ), activation dependent stiffness ( ), passive damping ( ), activation dependent damping ( ), visco-elasticity ( ) and activation dependent visco-elasticity ( ). This generalized model was fit to data from 7 different muscles, with amplitudes of motion ranging between 0.5 and 2.0 mm, and frequencies of oscillation ranging between 1 and 7 Hz. We found that this generalized model (Equation (S1)) captured 74.2%±5.5% of the temporal variance in the force generated.
To further simplify this model, and to assess the importance of each of those relative terms, we performed linear regression to fit the relationship between the observed force and the individual regressors ( ,˙ and ), and found that the dominant terms characterizing muscle response were the bilinear terms, namely , and . Figure S1C shows the contributions of each of the individual regressors to the variance of the output force. The results show that the bilinear approximation ≈ + + + captured 73.6% ± 5.7% of the total variance in muscle force, which is a minor degradation in accuracy compared to the model including the damping and visco-elastic terms. The parameters were normalized by the value of the coefficient to account for the differences in muscle size. The normalized parameters used in the optimization were: = 0, = 0.016, = 1, and = 0.727.
Fitting the parameters required knowledge of the activation . While estimates of from Equations (1) and (2) are possible, errors in the activation dynamics model may contribute to errors in the experimental identification of the impedance model. In order to minimize this effect, for each muscle we used a normalized profile of its isometric contraction as a representation of the activation dynamics. In doing so, we made the assumption that the activation dynamics did not vary with the length or velocity of the muscle. However, the normalized isometric force profile was used to determine the parameters of Equations (1) and (2). The parameters of the activation model used were = −62, = −62 and = 3947, which results in a step input gain of unity, and a twitch rise and fall time of 125 msec. We do not necessarily advocate that this model be used as a general muscle model since it has been tested only in the context of oscillatory motions with bursting contractions. However, in those particular conditions, it provided a reasonably descriptive approximation to the muscle force that is also mathematically tractable enough to be used for the optimal control formulation. What is evident from the model, however, is that the dependence of muscle force on velocity was found to be minimal. Whether this is a particular feature of that muscle, or whether it is a resultant effect of this particular experimental paradigm can be a topic of further research.

Computation of Optimal Control
The solution of the 2-point boundary value problem (Equation (7)) is non-trivial, primarily due to the discontinuous nature of the control. Even with the addition of the quadratic term, the hard bounds on the control lead to the failure of many of the standard solution algorithms. We solve this issue via the use of penalty methods. Penalty methods simplify the solution of a constrained optimization problem by converting it to a sequence of unconstrained optimization problems [2]. This sequence of problems is solved, with an increasing penalty function added for violating the constraints. In the limit, as the value of the penalty weight increases, the solutions will converge to the optimal solution of the constrained problem.
We applied penalty methods to the optimal control problem in the following manner. We would like to make sure that the control is limited to the admissible bounds, i.e.
≤ and ≥ for all . Therefore the penalties are imposed if the control violates those bounds. We also use quadratic penalties in the form of Therefore, the optimal control problem takes the following form: * = min Recall that the optimal control was a bang-bang control strategy, implying that depending on the sign of the switching function, it is optimal to push the control signal as far out as possible. With the penalty methods, a softened bound becomes optimal to push the control signal as far out until the penalty dominates.
Consider the case if the resulting optimal control is between and . Then ( ) = 0. Therefore, * 1 = arg min , then ( ) = ( − ) 2 . Therefore, * 1 = arg min In summary: * We solve this softened version of the problem numerically. We start with low values for , and increase it until ( ) violates the bounds constraints to within acceptable errors.
Solving the 2-point boundary value problem The 2-point boundary value problem is solved using the Matlab solver bvp4c for increasing values of . This solver descretizes the solution space to via an adaptive mesh, and fits piecewise cubic functions to each segment. The residuals in the interface points have to match the differential equation, and the boundary conditions. A description of the algorithm is found in [3].