Objectives

In the study of human movement, experimental measurement is generally limited to the kinematics of the body segments, external reaction forces, and electromyographic (EMG) signals. While these data are essential for characterizing movement, important information is missing. For example, because the body is actuated by more muscles than it has degrees of freedom, we cannot uniquely solve for the muscle forces that give rise to an observed motion. Yet, knowledge of muscle force is essential for quantifying the stresses placed on bones and also for understanding the functional roles of muscles in normal and pathological movement. Using dynamic models of the musculoskeletal system to simulate movement provides not only a means of estimating muscle forces, but also a framework for investigating how the various components of the musculoskeletal system interact to produce movement.

The purpose of this lab is to introduce you to the components of a musculoskeletal model, illustrate how these components can be combined, and demonstrate the value of dynamic simulation. You will use the results of dynamic simulations for feedback as you manually edit the excitation signals applied to the muscles of the lower extremity, with the goal of making a musculoskeletal model jump as high as possible. Jumping was chosen as the activity for this lab because it has a well-defined objective (i.e., jump as high as possible) and, although still complex, its muscular coordination is relatively simple compared to that of walking.

In this example, you will use several features of the OpenSim GUI to search for a set of muscle excitations that maximizes the jump height. In the course of this example, you will:

  • investigate the functions of muscles when they are fired in isolation and in conjunction with other muscle groups;
  • quantify the magnitude of the articular contact forces in the hip;
  • compare the ground contact forces produced by your simulation to published literature in the field; and
  • compare the force in muscles to the maximum isometric force during jumping.

Acknowledgements

This lab was originally created by Jeff Reinbolt, B.J. Fregly, Clay Anderson, Allison Arnold, Silvia Blemker, Darryl Thelen, and Scott Delp using torque actuators. Daniel Jacobs contributed to improving and refining the example.

Model

The Gait2392 and Gait2354 models are three-dimensional, 23-degree-of-freedom computer models of the human musculoskeletal system. The models were created by Darryl Thelen (University of Wisconsin-Madison) and Ajay Seth, Frank C. Anderson, and Scott L. Delp (Stanford University). The models feature lower extremity joint definitions adopted from Delp et al. (1990), low back joint and anthropometry adopted from Anderson and Pandy (1999), and a planar knee model adopted from Yamaguchi and Zajac (1989).

The Gait2392 model features 92 musculotendon actuators to represent 76 muscles in the lower extremities and torso. For the Gait2354 model, the number of muscles was reduced by Anderson to improve simulation speed for demonstrations and educational purposes. Seth removed the patella to avoid kinematic constraints; insertions of the quadriceps are handled with moving points in the tibia frame.

The default, unscaled version of these models represents a subject that is about 1.8 m tall and has a mass of 75.16 kg.

More details about how the models were constructed can be found on the Gait 2392 and 2354 Models page in the Musculoskeletal Models section.

Background

Dynamic models of the musculoskeletal system typically consist of four important components:

  1. The equations of motion for the body, or skeletal dynamics;
  2. A model of musculoskeletal geometry;
  3. A model of muscle–tendon mechanics; and
  4. A model of activation dynamics.



Figure 1. Forward Dynamic Simulation

Figure 1 illustrates how these components are combined to execute a forward dynamic simulation. Based on a set of initial states, which include the muscle activations \( \vec{a}\left(t\right)\), the muscle forces \( \vec{f}\left(t\right)\), the generalized speeds \( \dot{\vec{q}}\left(t\right)\), and the generalized coordinates \( \vec{q}\left(t\right)\), differential equations are used to compute the time rate of change of the states. At each time step of the simulation, the differential equations are numerically integrated to compute the states at the next time step. After each time step, the new states are used to calculate new derivatives and the forward dynamics process repeats, advancing the states in time until the final time of the simulation is reached. In the simulations you will conduct in this lab, a variable-step-size integrator is used.

Segment Dynamics

The equations of motion for the body allow one to compute the accelerations of the body segments when forces and torques are applied to the body. The equations of motion can be expressed as follows:

(1) \[ \ddot{\vec{q}} = I\left(\vec{q}\;\!\right)^{-1} \left( C\big(\vec{q}, \dot{\vec{q}}\,^2\big) + G\left(\vec{q}\;\!\right) + R\left(\vec{q}\;\!\right) \vec{f}_{M} + E\left(\vec{q}\;\!\right) \vec{f}_{E} \right)\]

Eq. (1) is simply an elaboration of Newton’s second law for a multi-link system, rearranged so that one can compute acceleration (i.e., \( \vec{a} = M\big(\vec{q}, \dot{\vec{q}}\;\!\big)^{-1} \vec{f}_M\)). The vector of generalized coordinates, \( \vec{q}\), is used to specify the position and orientation of the body segments. The time derivatives of \( \vec{q}\) (i.e., \( \dot{\vec{q}}\)and \( \ddot{\vec{q}}\)) represent the velocities and accelerations of the segments.

Depending on how one chooses to model the body, elements of \( \vec{q}\) may be translations, orientations of segments with respect to the lab frame (segment angles), or orientations of segments with respect to other segments (joint angles). Implicit in one’s choice of generalized coordinates are one’s assumptions about how the joints of the body function. For example, one often models the hip joint as a three-degree-of-freedom ball-and-socket (spherical) joint, which requires three generalized coordinates: flexion-extension (\( q_1\)), abduction-adduction (\( q_2\)), and internal-external rotation (\( q_3\)).

The system mass matrix, \( I\left(\vec{q}\;\!\right)\), characterizes the inertial properties of the body (i.e., masses and moments of inertia). The remaining terms in Eq. (1) express the generalized forces or torques that act on the body segments. \( C\big(\vec{q}, \dot{\vec{q}}\,^2\big)\) represents centripetal forces that arise from the angular velocities of the segments, \( G\left(\vec{q}\;\!\right)\) represents gravitational forces, \( R\left(\vec{q}\;\!\right) \vec{f}_{M}\) represents the moments applied at the joints by the muscles, and \( E\left(\vec{q}\;\!\right) \vec{f}_{E}\) represents external forces applied to the body, such as the ground reaction force. The matrix \( R\left(\vec{q}\;\!\right)\) is a matrix of moment arms that transform the muscle forces, \( \vec{f}_{M}\), into joint torques. The matrix \( E\left(\vec{q}\;\!\right)\) performs a similar function for the external forces, \( \vec{f}_{E}\).

For simple models, it is possible to derive the equations of motion by hand; however, for more complex models, this is generally not feasible (and it is highly error-prone!), and the equations of motion are generated on a computer. The jumping model used in this lab has 23 degrees of freedom (Anderson and Pandy, 1999), and the equations of motion for the jumping model were generated using OpenSim (Delp et al., 2007).

 

Musculoskeletal Geometry

Accurately representing the path of a muscle from its origin to its insertion is one of the more challenging aspects of modeling the musculoskeletal system. Sometimes a muscle can be represented as a straight-line path between its origin and insertion. Other times, it is adequate to approximate the path as a series of straight line segments that pass through a series of via points (Delp et al., 2007). When modeling muscle paths in three dimensions, it is often necessary to simulate how muscles wrap over underlying bone or musculature. Cylinders, spheres, and ellipsoids have been used as wrapping surfaces (Van der Helm et al., 1992; Garner and Pandy, 2000; Arnold et al., 2000) (see figure at left).

 

 

 


Muscle–Tendon Mechanics

The force-producing properties of muscle are complex and nonlinear (see McMahon (1984) for review) (Fig. 1). For simplicity, lumped-parameter, dimensionless muscle models, capable of representing a range of muscles with different architectures, are most commonly used in the dynamic simulation of movement (Zajac, 1989). In a complex musculoskeletal model, the model can be actuated by 50 or more muscle–tendon units, each of which is represented as a Hill-type contractile element in series with a tendon.

The Thelen2003Muscle model uses a standard equilibrium muscle model based on the Hill model. The muscle–tendon complex consists of three components: a contractile element (CE), a parallel element (PE), and a series element (SE). The muscle force generated is a function of three factors: the activation value (a), the normalized length of the muscle unit, and the normalized velocity of the muscle unit. The functions describing the force generated by a muscle as its length varies are called the active length curve (AL) for the contractile element and the passive length curve (PL) for the parallel element.

The parameters used to characterize each muscle are maximum isometric force, optimal muscle fiber length, tendon slack length, maximum contraction velocity, and pennation angle. For a table of muscle and tendon parameters, consult Anderson and Pandy (1999). During a forward dynamic simulation, the muscle force is calculated using two states: the activation value and the muscle fiber length.

 

Figure 1 from Thelen (2003).

 

Activation Dynamics

A muscle is not capable of generating force or relaxing instantaneously. The development of force is a complex sequence of events that begins with the firing of motor units and culminates in the formation of actin–myosin cross-bridges within the myofibrils of the muscle. When the motor units of a muscle depolarize, action potentials are elicited in the fibers of the muscle and cause calcium ions to be released from the sarcoplasmic reticulum. The increase in calcium ion concentrations then initiates the cross-bridge formation between the actin and myosin filaments (see Guyton (1986) for review). In isolated muscle twitch experiments, the delay between a motor unit action potential and the development of peak force has been observed to vary from as little as 5 milliseconds for fast ocular muscles to as much as 40 or 50 milliseconds for muscles comprised of higher percentages of slow-twitch fibers. The relaxation of muscle depends on the re-uptake of calcium ions into the sarcoplasmic reticulum. This re-uptake is a slower process than the calcium ion release, and so the time required for muscle force to fall can be considerably longer than the time for it to develop.

In the forward dynamic simulations you will conduct in this lab, activation dynamics is modeled using a first-order differential equation to relate the rate of change in activation (i.e., the concentration of calcium ions within the muscle) to excitation (i.e., the firing of motor units):

(2) \[ \dot{a} = \frac{x^2-xa}{\tau_\textrm{rise}} + \frac{x-a}{\tau_\textrm{fall}}\]

where \( a\) is the activation level of a muscle, \( x\) is the excitation level of a muscle, and \( \tau_\textrm{rise}\) and \( \tau_\textrm{fall}\) are the rise and fall time constants for activation, respectively. In the model, activation is allowed to vary continuously between zero (no contraction) and one (full contraction). In the body, the excitation level of a muscle is a function of both the number of motor units recruited and the firing frequency of the motor units. Some models for excitation–contraction coupling distinguish between these two control mechanisms (Hatze, 1976), but it is often not computationally feasible to use such models when conducting complex dynamic simulations.

Ground Contact Dynamics

The contact dynamics between the model bodies and the ground are represented with a continuous, nonlinear contact model published by Hunt and Crossley (1975). The contact force is modeled as

(3) \[ f_\textrm{hc} = - \lambda\dot{x}x^n - kx^n\]

where \( x\) and \( \dot{x}\) are the interference and interference rate between the contacting surfaces, \( \lambda\) is the damping constant, \( k\) is the spring constant, and \( n\) is the power constant. The contact surface between the foot and the ground is modeled as five spheres attached to the calcaneus and toe bodies on each foot.

Getting Started

A copy of these instructions can be saved as a .pdf or a Word document by selecting "Export to PDF" or "Export to Word" from the Tools menu in the top-right corner of this page. All the necessary files for this example are attached to the instructional page Coordinating Muscles for Optimal Jump Performance.

If you are completing this example as a laboratory exercise for a course on human movement, you will need to submit answers to the questions on the Questions: Optimal Jump Performance page.

References

  1. Anderson, F.C., Pandy, M.G. (1999). A dynamic optimization solution for vertical jumping in three dimensions. Computer Methods in Biomechanics and Biomedical Engineering, 2(3):201–231.
  2. Arnold, A.S., Salinas, S., Asakawa, D.J., Delp, S.L. (2000). Accuracy of muscle moment arms estimated from MRI-based musculoskeletal models of the lower extremity. Computer Aided Surgery, 5(2):108–119.
  3. Atkinson, L.V., Harley, P.J., Hudson, J.D. (1989). Numerical Methods with FORTRAN 77: A Practical Introduction. Addison–Wesley Publishing Company, Menlo Park.
  4. Delp, S.L., Loan, J.P., Hoy, M.G., Zajac, F.E., Topp, E.L., Rosen, J.M. (1990). An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures. IEEE Transactions on Biomedical Engineering, 37(8):757–767.
  5. Garner, B.A., Pandy, M.G. (2000). The obstacle-set method for representing muscle paths in musculoskeletal models. Computer Methods in Biomechanics and Biomedical Engineering, 3(1):1–30.
  6. Guyton, A.C. (1986). Textbook of Medical Physiology, 7th ed. W. B. Saunders Company, Philadelphia.
  7. Hatze, H. (1976). The complete optimization of a human motion. Mathematical Biosciences, 28(1–2):99–135.
  8. McMahon, T.A. (1984). Muscles, Reflexes, and Locomotion. Princeton University Press, Princeton.
  9. Symbolic Dynamics, Inc. (1996). SD/FAST User’s Manual, Version B.2. Mountain View.
  10. Van der Helm, F.C.T., Veeger, H.E.J., Pronk, G.M., Van der Woude, L.H.V., Rozendal, R.H. (1992). Geometry parameters for musculoskeletal modeling of the shoulder system. Journal of Biomechanics, 2:129–144.
  11. Zajac, F.E. (1989). Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. CRC Critical Reviews in Biomedical Engineering (Edited by JR Bourne), 17(4):359–411.
  12. Hunt, K.H., Crossley, F.R.E. (1975). Coefficient of restitution interpreted as damping in vibroimpact. Journal of Applied Mechanics, 42(2):440–445.