Muscle Optimizer Theory

Overview

Musculoskeletal modeling is a computational technique that allows simulating human motion and estimating internal joint loads and muscle forces. The aim is to gain deeper insight on the forces occurring in the human body. Simulations can be created using generic models derived from cadaveric data or using personalized models of the musculoskeletal geometry built from medical images.

However, these personalized models are currently challenging to generate because, while personalized bone geometries and joint parameters can be obtained from segmentation of medical images, musculotendon parameters that regulate force generation currently cannot be easily measured or estimated. This is especially true for the optimal fiber length \(l_o^m\) and tendon slack length \(l_s^t\).

The algorithms available in the literature to estimate these parameters can be divided into “anthropometric” and “functional” approaches, depending if they use as input only skeletal dimensions or integrate those with additional measurements of muscular function.

The anthropometric algorithm implemented in the Muscle Optimizer toolbox maps the muscle fiber operating ranges of a generic model into a scaled model, similarly to [2], but using a generic and efficient implementation. Musculoskeletal models to be optimized, computational muscle models and involved degrees of freedom are fully customizable by the user.

Theory

Dimensionless muscle model

The Muscle Optimizer toolbox works on OpenSim [3] models including dimensionless Hill-type muscle actuators of the kind proposed by Zajac in [4] (Fig. 1). These models include a contractile element (CE) connected to an elastic element in series (SE) and one in parallel (PE) and use normalized functions to define the active and passive muscle force-length curves, force-velocity curve and tendon force-strain curve. They need five parameters to be defined: optimal fiber length \(l_o^m\), tendon slack length \(l_s^t\), maximum fiber isometric force \(F_{iso}\)., pennation angle at optimal fiber length \(\alpha_o\) and maximum contraction velocity \(v_{max}\). Equilibrium of the musculotendon unit is obtained when the tendon and the fiber forces along the tendon direction are equal.

Fig. 1 Hill muscle model: A) Representation of the three elements Hill-type muscle model assumed by the software. B) Generic curves defining the dimensionless material properties of the tendon (left side, dashed line identifying tendon stretch for which tensile force equals maximum isometric force) and muscle (right side, solid line: active force, dashed line: passive force, dash-dot line: total muscle force). An equilibrated contraction of a non-pennate muscle is also represented (red dashed line).

Describing the dimensionless muscle operating range

The entire operating range and isometric force generating modalities of a muscle can be described by its normalized coefficients calculated as a function of the joint angles spanned by the muscle in the model. In this toolbox, we tackled this problem following the procedure described in our recent paper [1] and reported below.

The length of a musculotendon actuator \(l^{mt}\) can be calculated from the muscle length \(l^m\) and tendon length \(l^t\) as follows (Fig. 1 A):

(1)\[l^{mt} = l^m \cos\alpha + l^t\]

where \(\alpha\) is the pennation angle at that specific muscle length, calculated assuming constant muscle thickness, as:

(2)\[\alpha = \arcsin(\frac{l_o^m \sin \alpha_o}{l^m})\]

Following [4] we defined the normalized fiber length as:

(3)\[\tilde{l}^m = \frac{l^m}{l^m_o}\]

and we decided to define a normalized tendon length (\(\epsilon^t\) being tendon strain) as:

(4)\[\tilde{l}^t = \frac{l^t}{l^t_s}=(1+\epsilon^t)\]

Using normalized coefficients, the musculotendon length can now be expressed as \(f(l_o^m, l_s^t)\):

(5)\[l^{mt} = (\tilde{l}^m \cos\alpha)l_o^m + \tilde{l}^t l^t_s\]

If the normalized coefficients are calculated from \(l^m\) and \(l^t\) in a musculotendon unit equilibrated for isometric contraction at maximum activation, Eq. (5) ensures static equilibrium between the tendon and (active plus passive) muscle force of any muscle actuator of length \(l^{mt}\) whose \(l_o^m\) and satisfy that equation. For an example, see Fig. 1 B.

Musculotendon parameters estimation

Fig. 2 Flowchart from [1] representing the main steps of the proposed algorithm for estimating the parameters of each musculotendon unit (blue: reference model, red: target model, green: optimized model).

Similar to [2], the proposed method aims to map the normalized muscle operating conditions of an existing “reference model”, whose muscle parameters are assumed to be physiologically valid, onto a “target model” of different anthropometric dimensions for equivalent joint configurations. The algorithm, represented as flow chart in Fig. 2 A, consists of the following steps, applied to each muscle included in the model:

1. In the reference model, the \(N_q\) degrees of freedom spanned by the musculotendon actuator are uniformly sampled using \(n_{dof}\) points per coordinate, so generating a set of \(n=(n_{dof})^{N_q}\) total model poses in which the considered muscle is equilibrated for isometric contraction with a unitary activation.
2. For each pose \(i=1,2,...n\) of the reference model, the vectors of pennation angles \(\alpha_{i,ref}\), normalized fiber lengths \(\tilde{l}^m_{i,ref}\) and normalized tendon lengths \(\tilde{l}^t_{i,ref}\) are calculated using Eqs. (2) (3) (4). Normalized fiber lengths outside the range of 0.5 and 1.5 or causing pennation angles approaching 90 degrees are excluded.
3. The musculotendon lengths \(\alpha_{i,ref}\) are calculated in the target model for the same joint configurations used in the reference model at step 1.
4. Imposing Eq. (5) for all poses, the following \(n \times 2\) linear system is obtained and can be solved in a least square sense, so obtaining the unknown values of \(l^m_{o,targ}\) and \(l^t_{s,targ}\).

(6)\[\begin{split}\begin{bmatrix} l_1^{mt} \\ l_2^{mt} \\ \vdots \\ l_n^{mt} \end{bmatrix}_{targ} = \begin{bmatrix} \tilde{l}_1^m\cos{\alpha_1} & \tilde{l}_1^t \\ \tilde{l}_2^m\cos{\alpha_2} & \tilde{l}_2^t \\ \vdots \\ \tilde{l}_n^m\cos{\alpha_n} & \tilde{l}_n^t \end{bmatrix}_{ref} \begin{bmatrix} l_o^m \\ l_s^t \end{bmatrix}_{targ}\end{split}\]

5. If the system yields non physiological results, such as \(l^t_{s,targ}<0\), the tendon fraction of \(l^mt_{targ}\) is temporary fixed at the same proportion of the reference model \(l^m_{o,targ}\) estimated by solving Eq. (6) and finally \(l^t_{s,targ}\) is re-computed from the same equation using the new value of \(l^m_{o,targ}\).

Bibliography

[1]	(1, 2) Luca Modenese, Elena Ceseracciu, Monica Reggiani, and David G Lloyd. Estimation of musculotendon parameters for scaled and subject specific musculoskeletal models using an optimization technique. Journal of biomechanics, 49(2):141–148, 2016.

[2]	(1, 2) C.R. Winby, D.G. Lloyd, and T.B. Kirk. Evaluation of different analytical methods for subject-specific scaling of musculotendon parameters. Journal of Biomechanics, 41(8):1682 – 1688, 2008.

[3]	S. L. Delp, F. C. Anderson, A. S. Arnold, P. Loan, A. Habib, C. T. John, E. Guendelman, and D. G. Thelen. Opensim: open-source software to create and analyze dynamic simulations of movement. IEEE Transactions on Biomedical Engineering, 54(11):1940–1950, Nov 2007.

[4]	(1, 2) Felix E Zajac. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical reviews in biomedical engineering, 17(4):359–411, 1989.