|
Body segments are connected to other segments at joints, which represent articulations (e.g. hip, knee, ankle) within the body. Joints are frictionless and non-deformable and, consequently, forces and moments transfer across these joints in an equal and opposite manner. A joint definition includes its location relative to the body segments it connects, as well as the degrees of freedom (DOFs) permissible at the joint (one, two or three axes of rotations and, less commonly, translations). The direction of rotation axes are usually expressed in terms of a local (body segment) coordinate system.
신체 세그먼트는 신체 내의 관절(예: 엉덩이, 무릎, 발목)을 나타내는 관절에서 다른 세그먼트와 연결됩니다. 관절은 마찰이 없고 변형이 불가능하며, 따라서 힘과 모멘트는 관절을 통해 동등하고 반대되는 방식으로 전달됩니다. 관절의 정의에는 관절이 연결된 신체 세그먼트에 대한 위치와 관절에서 허용되는 자유도(DOF)(1, 2 또는 3축의 회전축 및 덜 일반적으로는 이동축)가 포함됩니다. 회전축의 방향은 일반적으로 로컬(바디 세그먼트) 좌표계로 표현됩니다.
With the body segments and joints defined, the skeletal portion of the musculoskeletal model is sufficient to conduct dynamic analyses. Musculoskeletal models extend the analysis by including models of individual muscles or even portions of muscles (figure 1). Muscle models range in their complexity and, as with other parts of the model, the necessary complexity depends on the motion being simulated. (For a description of how muscle complexity can be modelled, see §5.3.) A muscle model minimally requires geometry (i.e. path definition) and a parameter representing its maximum possible force. More complex muscle models also include parameters for pennation angle, optimal fibre length, physiological cross-sectional area (PCSA) and activation and deactivation times. Muscle geometry must include the coordinates of its origin (location relative to the proximal segment) and insertion (location relative to the distal segment), but muscles can be geometrically more complex than a straight line path between origin and insertion. The most basic of these complexities derives from a muscle with a broad origin (e.g. gluteus medius) or insertion (e.g. adductor magnus). In these muscles, a single origin-to-insertion path does not encompass their possible actions. One way to solve this is to split the muscle into muscle parts, referred to here as muscle elements. Each muscle element represents a region of the muscle (e.g. the superior portion of adductor magnus) and can have its own set of muscle parameters.
Most modelling systems allow for additional coordinate positions that are used to define the curved path that the physical muscles take between bony attachments (e.g. the path sartorius takes across the anterior thigh as it traverses from anterior superior iliac spine (ASIS) to proximal tibia). Via points, the simplest type of path point, are one or more coordinate positions that are intermediate between the origins and insertion points. The muscle assumes a series of straight line paths between the origin, via point(s) and insertion. The location of the via point is fixed relative to the segment's centre of mass and the connection to other segments. With via points, the muscle is required to pass through each specific point at all times. These points allow muscles to pass over other body surfaces (e.g. iliacus wrapping over the pelvic brim). In addition, muscle geometry paths can be altered by the use of ‘wrapping surfaces’ that are designed to account for relative motion between musculoskeletal elements and other soft tissue (e.g. rectus femoris sliding across the soft tissue structures anterior to the hip). These wrapping surfaces are geometric objects (e.g. cylinders, ellipsoids) that are fixed relative to a particular body segment. Muscle paths are not constrained to pass through any particular coordinate point (as with via points), but rather can move freely over the wrapping surface to find the shortest possible path across the surface between origin, insertion or via points that are beyond the wrapping surface. As the proximal and distal segments change orientation relative to each other, the point of contact between the muscle and the adjacent musculoskeletal elements can move (e.g. gastrocnemius relative to the posterodistal femur as the knee flexes). Forces developed from the action of the muscle act on the body segment used to define the geometry (origin, insertion, via points and/or points on wrapping surfaces). Muscles may also be modelled to be more complex in terms of the number of parameters that define force generation. Minimally, a muscle model must include a maximum force value which defines the force when the muscle is maximally activated, but other approximations of muscle force generation are available, e.g. the Hill-type muscle model [64].
The final required model element is a set of virtual motion markers, which is needed to match model motion to motion measured experimentally (figure 1). For motion tracking (described below), the human subject is fitted with a set of markers that are tracked by a camera system. The musculoskeletal model must include one virtual marker for each experimentally measured marker. Each virtual marker is connected to a body segment and defined in terms of that segment's coordinate system. For instance, if researchers affix a motion-tracking marker over the ASIS of a human participant, the model will require a virtual marker connected to the pelvis body segment with coordinates that position it correctly relative to the model ASIS.
신체 세그먼트와 관절이 정의되면
근골격 모델의 골격 부분만으로도
동적 분석을 수행하기에 충분합니다.
근골격 모델은
개별 근육 또는 근육의 일부 모델을 포함함으로써 분석을 확장합니다(그림 1).
근육 모델은 복잡도가 다양하며,
모델의 다른 부분과 마찬가지로 필요한 복잡도는 시뮬레이션되는 동작에 따라 달라집니다.
(근육 복잡도를 모델링하는 방법에 대한 설명은 §5.3을 참조하세요.)
근육 모델에는 최
소한의 지오메트리(즉, 경로 정의)와 최대 가능한 힘을 나타내는
파라미터가 필요합니다.
더 복잡한 근육 모델에는
페네이션 각도,
최적의 근섬유 길이,
생리적 단면적(PCSA),
활성화 및 비활성화 시간에 대한 파라미터도 포함됩니다.
More complex muscle models also include parameters for
pennation angle,
optimal fibre length,
physiological cross-sectional area (PCSA) and
activation and deactivation times.
근육 지오메트리에는 원점(근위 세그먼트 기준 위치)과 삽입(원위 세그먼트 기준 위치)의 좌표가 포함되어야 하지만, 근육은 원점과 삽입 사이의 직선 경로보다 기하학적으로 더 복잡할 수 있습니다. 이러한 복잡성 중 가장 기본적인 것은 기원이 넓은 근육(예: 중둔근) 또는 삽입부(예: 내전근)에서 비롯됩니다. 이러한 근육에서는 단일 기원-삽입 경로가 가능한 동작을 모두 포괄하지 못합니다. 이 문제를 해결하는 한 가지 방법은 근육을 근육 요소라고 하는 근육 부분으로 분할하는 것입니다. 각 근육 요소는 근육의 영역(예: 내전근의 상부 부분)을 나타내며 자체적인 근육 파라미터 세트를 가질 수 있습니다.
대부분의 모델링 시스템에서는 실제 근육이 뼈 부착물 사이에서 곡선 경로를 정의하는 데 사용되는 추가 좌표 위치를 허용합니다(예: 전방 장골(ASIS)에서 근위 경골로 이동할 때 허벅지 앞쪽을 가로지르는 사르토리우스 근육의 경로). 가장 단순한 유형의 경로 지점인 비아 포인트는 원점과 삽입점 사이의 중간 지점인 하나 이상의 좌표 위치입니다. 근육은 원점, 경유점, 삽입점 사이의 일련의 직선 경로를 가정합니다. 비아 포인트의 위치는 세그먼트의 질량 중심과 다른 세그먼트와의 연결을 기준으로 고정됩니다. 비아 포인트를 사용하면 근육은 항상 각 특정 지점을 통과해야 합니다. 이러한 지점을 통해 근육이 다른 신체 표면을 통과할 수 있습니다(예: 장요근이 골반 가장자리를 감싸고 있는 경우). 또한 근골격 요소와 다른 연조직 사이의 상대적인 움직임을 고려하도록 설계된 '랩핑 표면'을 사용하여 근육 지오메트리 경로를 변경할 수 있습니다(예: 엉덩이 앞쪽의 연조직 구조를 가로질러 미끄러지는 대퇴직근). 이러한 래핑 표면은 특정 신체 세그먼트에 대해 고정된 기하학적 개체(예: 원통, 타원체)입니다. 근육 경로는 비아 포인트와 같이 특정 좌표점을 통과하도록 제한되지 않고 랩핑 표면 위에서 자유롭게 이동하여 원점, 삽입점 또는 랩핑 표면 너머의 비아 포인트 사이에서 표면을 가로지르는 가능한 최단 경로를 찾을 수 있습니다. 근위 및 원위 세그먼트가 서로에 대해 방향을 변경함에 따라 근육과 인접한 근골격 요소 사이의 접촉 지점이 움직일 수 있습니다(예: 무릎이 구부러질 때 대퇴골 후방에 대한 비복근). 근육의 작용으로 인해 발생하는 힘은 지오메트리를 정의하는 데 사용되는 신체 세그먼트에 작용합니다(원점, 삽입, 경유점 및/또는 래핑 표면의 포인트). 근육은 힘 생성을 정의하는 파라미터의 수 측면에서 더 복잡하게 모델링할 수도 있습니다. 근육 모델에는 최소한 근육이 최대로 활성화되었을 때의 힘을 정의하는 최대 힘 값이 포함되어야 하지만, 언덕형 근육 모델[64]과 같이 근육 힘 생성에 대한 다른 근사치를 사용할 수 있습니다.
마지막으로 필요한 모델 요소는 가상 모션 마커 세트이며, 이는 모델 모션을 실험적으로 측정된 모션과 일치시키는 데 필요합니다(그림 1). 모션 추적(아래 설명)을 위해 사람 피사체에는 카메라 시스템으로 추적되는 마커 세트가 장착됩니다. 근골격계 모델에는 실험적으로 측정된 각 마커에 대해 하나의 가상 마커가 포함되어야 합니다. 각 가상 마커는 신체 세그먼트에 연결되며 해당 세그먼트의 좌표계에 따라 정의됩니다. 예를 들어 연구자가 인간 참가자의 ASIS 위에 동작 추적 마커를 부착하는 경우, 모델에는 골반 신체 세그먼트에 연결된 가상 마커가 모델 ASIS를 기준으로 정확한 위치를 지정하는 좌표와 함께 필요합니다.
2.2. The input data: motion data and ground reaction forces
An inverse dynamics approach to musculoskeletal modelling (described below) requires that the motions (translation and rotation) of the rigid segments (i.e. the skeletal elements) and externally applied forces be known. This process is well established and has been described in detail elsewhere (e.g. [65–67]). To obtain motions, retroreflective (infrared) markers attached to the skin or clothing are frequently used. Each body segment requires a minimum of three markers that are trackable throughout the entire sequence of movements, so this traceability should be checked in all portions of the data collection volume. The three-dimensional ground reaction force (GRF) for a complete stance is also necessary. GRFs are typically acquired through force plates embedded either in the floor or in a treadmill.
2.2. 입력 데이터: 모션 데이터 및 지면 반력
근골격 모델링에 대한 역역학 접근법(아래 설명)을 사용하려면 강체 세그먼트(즉, 골격 요소)의 움직임(이동 및 회전)과 외부에서 가해지는 힘을 알아야 합니다. 이 프로세스는 잘 확립되어 있으며 다른 곳에서 자세히 설명되어 있습니다(예: [65-67]). 동작을 얻기 위해 피부나 옷에 부착된 역반사(적외선) 마커가 자주 사용됩니다. 각 신체 세그먼트에는 전체 움직임 시퀀스에서 추적 가능한 최소 3개의 마커가 필요하므로 데이터 수집 볼륨의 모든 부분에서 이 추적 가능성을 확인해야 합니다. 완전한 자세를 위한 3차원 지면 반작용력(GRF)도 필요합니다. GRF는 일반적으로 바닥이나 러닝머신에 내장된 포스 플레이트를 통해 수집합니다.
2.3. The simulation
2.3.1. Scaling and marker optimization
Model scaling and virtual marker position optimization are the processes of determining the model segment dimensions and virtual marker positions so that the virtual model matches the human subject and experimental protocol (motion marker placement on the human subject) as well as possible [68]. The degree of matching between the model and experimental data is evaluated using a global error metric between experimental and virtual motion marker positions. During the inverse kinematic analysis (described below), musculoskeletal modelling software finds the position of the model, including overall location in space and individual joint positions, that minimizes the sum of squared deviations between the virtual marker set and experimental motion marker data. The goal of marker optimization and body scaling is to minimize this error term.
The human models that are included with many modelling software packages are generic human models. For analysis, the model must match the subject for whom the input data are known in terms of mass, limb parameters (e.g. distances between joint centres and segment moments of inertia) and virtual marker positions [69,70]. Marker position optimization and scaling are intertwined processes because marker positions can, in part, determine segment dimensions. Initial scaling of the generic virtual model achieves a loose match between the overall bodily dimensions of the subject and those of the virtual model. If a virtual model from a repository is used (rather than creating a virtual model from scratch), then the model will have been developed using the parameters of the individual whose data were available (i.e. a generic model). This generic model will almost certainly not represent the dimensions of the subject on whom motion-tracking data and GRFs were collected, so the generic model needs to be scaled [69]. The simplest scaling approach is uniform scaling of lengths and volumes. Total body mass is the most basic dimension needed, and the ratio of the body mass of the subject to that of the generic model is typically used to adjust segment masses and other parameters that are volume-dependent. The ratio of subject to generic model stature, a length scale, is used as an initial estimation of segment dimensions and other dimensions that are length-dependent. Other anthropometric data (e.g. distances between landmarks) can be used instead of, or in addition to, stature. Other scaling methods are also available, including fully manual scaling, using body fat to adjust mass distribution among the segments, and hybrid scaling (e.g. manually scaling some parameters and algorithmically scaling others). Furthermore, depending on the sophistication of the software used to complete the mechanical analysis, these initial scalings may be adjusted during marker optimization.
Different software packages conduct marker optimization and secondary scaling of segments using different algorithms and require varying levels of interaction/input from the user. At one end of the user-input spectrum, some systems require that the user optimizes the locations of the markers on the model ‘by hand’, moving them relative to their parent segment and then interrogating error values associated with each marker during motion tracking. This process is tedious and may require significant effort before optimal marker position is achieved. At the other end of the spectrum, some software packages change both segment dimensions and virtual motion-tracking markers automatically in a single cohesive optimization process that requires little user interaction. In these frameworks, the user can permit or restrict changes in virtual marker position in none, some or all of the coordinate directions. For instance, motion-tracking markers associated with easily palpable bony landmarks (e.g. medial malleolus) are also easy to locate on the virtual model (referencing the virtual skeleton) and may be considered ‘fixed’ in their location relative to their body segment (the shank). Thus, during the optimization process, such a fixed marker would not change its position relative to the segment, but would influence the determination of shank length. Other markers (e.g. marker plates attached mid-segment) may be allowed to optimize their position in all three coordinate directions, and then contribute little to segment length determination. Finally, some markers may be well defined in two coordinate directions, but less so in the third. For instance, the superoinferior and mediolateral positions of an ASIS marker may be well defined and constrained from moving relative to the pelvis during the optimization process. The anteroposterior position of the virtual ASIS marker relative to the virtual pelvis may be more difficult to estimate as it reflects tissue thickness over this bony landmark (which is hard to measure on living subjects in the absence of internal imaging modalities). Thus, the user might fix the ASIS marker in two coordinate directions, but allow its position to optimize in the third. The fixed mediolateral and superoinferior position on the two ASISs would then influence the scaling of pelvic width.
2.3. 시뮬레이션
2.3.1. 스케일링 및 마커 최적화
모델 스케일링 및 가상 마커 위치 최적화는 모델 세그먼트 치수와 가상 마커 위치를 결정하여 가상 모델이 인체 피험자 및 실험 프로토콜(인체 피험자에 대한 모션 마커 배치)과 최대한 일치하도록 하는 프로세스입니다[68]. 모델과 실험 데이터 간의 일치 정도는 실험과 가상 모션 마커 위치 사이의 글로벌 오차 지표를 사용하여 평가합니다. 역 운동학 분석(아래 설명)을 수행하는 동안 근골격계 모델링 소프트웨어는 공간에서의 전체 위치와 개별 관절 위치를 포함하여 가상 마커 세트와 실험 모션 마커 데이터 간의 제곱 편차 합을 최소화하는 모델의 위치를 찾습니다. 마커 최적화 및 바디 스케일링의 목표는 이 오차 항을 최소화하는 것입니다.
많은 모델링 소프트웨어 패키지에 포함된 인간 모델은 일반적인 인간 모델입니다. 분석을 위해 모델은 질량, 사지 매개변수(예: 관절 중심과 세그먼트 관성 모멘트 사이의 거리) 및 가상 마커 위치 측면에서 입력 데이터가 알려진 피험자와 일치해야 합니다[69,70]. 마커 위치가 부분적으로 세그먼트 치수를 결정할 수 있기 때문에 마커 위치 최적화와 스케일링은 서로 얽혀 있는 프로세스입니다. 일반 가상 모델의 초기 스케일링은 피험자의 전체 신체 치수와 가상 모델의 치수가 느슨하게 일치하도록 합니다. 처음부터 가상 모델을 생성하지 않고 저장소의 가상 모델을 사용하는 경우, 해당 모델은 데이터를 사용할 수 있는 개인의 매개변수(즉, 일반 모델)를 사용하여 개발되었을 것입니다. 이 일반 모델은 동작 추적 데이터와 GRF가 수집된 피험자의 차원을 거의 정확하게 나타내지 못하므로 일반 모델의 크기를 조정해야 합니다[69]. 가장 간단한 스케일링 접근 방식은 길이와 부피를 균일하게 스케일링하는 것입니다. 총체질량은 필요한 가장 기본적인 차원이며, 일반적으로 피험자의 체질량과 일반 모델의 체질량 비율을 사용하여 체적에 따라 달라지는 세그먼트 질량 및 기타 파라미터를 조정합니다. 피험자 대 일반 모델 키의 비율, 즉 길이 척도는 세그먼트 치수 및 길이에 따라 달라지는 기타 치수의 초기 추정치로 사용됩니다. 키 대신 또는 키와 함께 다른 인체 측정 데이터(예: 랜드마크 사이의 거리)를 사용할 수 있습니다. 완전 수동 스케일링, 체지방을 사용하여 세그먼트 간의 질량 분포 조정, 하이브리드 스케일링(예: 일부 매개변수는 수동으로 스케일링하고 다른 매개변수는 알고리즘으로 스케일링) 등 다른 스케일링 방법도 사용할 수 있습니다. 또한 기계적 분석을 완료하는 데 사용되는 소프트웨어의 정교함에 따라 이러한 초기 스케일링은 마커 최적화 중에 조정될 수 있습니다.
소프트웨어 패키지마다 서로 다른 알고리즘을 사용하여 마커 최적화와 세그먼트의 2차 스케일링을 수행하며 사용자의 다양한 수준의 상호 작용/입력을 요구합니다. 사용자 입력 스펙트럼의 한쪽 끝에서, 일부 시스템은 사용자가 모델에서 마커의 위치를 '손으로' 최적화하여 부모 세그먼트를 기준으로 이동한 다음 모션 추적 중에 각 마커와 관련된 오류 값을 조사하도록 요구합니다. 이 프로세스는 지루하고 최적의 마커 위치를 얻기까지 상당한 노력이 필요할 수 있습니다. 반면에 일부 소프트웨어 패키지는 사용자 상호작용이 거의 필요 없는 단일 통합 최적화 프로세스에서 세그먼트 치수와 가상 모션 트래킹 마커를 모두 자동으로 변경합니다. 이러한 프레임워크에서는 사용자가 가상 마커 위치의 좌표 방향 중 일부 또는 전부에 대한 변경을 허용하거나 제한할 수 있습니다. 예를 들어, 쉽게 만져지는 뼈의 랜드마크(예: 내측 연골)와 관련된 동작 추적 마커는 가상 모델에서 쉽게 찾을 수 있으며(가상 골격을 참조), 해당 신체 세그먼트(정강이)를 기준으로 그 위치가 '고정'된 것으로 간주될 수 있습니다. 따라서 최적화 프로세스 중에 이러한 고정 마커는 세그먼트에 대한 위치가 변경되지 않지만 생크 길이 결정에 영향을 미칩니다. 다른 마커(예: 세그먼트 중간에 부착된 마커 플레이트)는 세 좌표 방향 모두에서 위치를 최적화할 수 있지만 세그먼트 길이 결정에 거의 기여하지 않을 수 있습니다. 마지막으로, 일부 마커는 두 좌표 방향에서는 잘 정의되지만 세 번째 좌표 방향에서는 그렇지 않을 수 있습니다. 예를 들어, ASIS 마커의 상하 및 중측 위치는 잘 정의되어 있지만 최적화 프로세스 중에 골반을 기준으로 움직이지 못하도록 제한될 수 있습니다. 가상 골반에 대한 가상 ASIS 마커의 전후방 위치는 이 뼈 랜드마크 위의 조직 두께를 반영하기 때문에 추정하기가 더 어려울 수 있습니다(내부 이미징 장비가 없는 경우 살아있는 피사체에서 측정하기 어려움). 따라서 사용자는 ASIS 마커를 두 좌표 방향으로 고정하되 세 번째 좌표 방향에서 위치가 최적화되도록 할 수 있습니다. 그러면 두 ASIS의 고정된 내측 및 상측 위치가 골반 폭의 스케일링에 영향을 미칩니다.
2.3.2. Inverse kinematics/motion tracking
To carry out inverse dynamics and, in turn, estimate internal forces and moments, segment linear and rotational accelerations must be obtained. Most musculoskeletal modelling software (e.g. Anybody Modeling System™, OpenSim) uses a kinematic analysis procedure of over-determined biomechanical systems, which is formulated as a weighted least-squares optimization problem that matches experimental and model marker positions (e.g. [47,71]).
2.3.3. Inverse dynamics
While statics is the branch of mechanics that deals with stationary objects, dynamics is the branch in which motions are the focus. In statics, an equilibrium that greatly facilitates examination of the system exists where the environment (e.g. the substrate) and the structure (e.g. the human body) are balanced. In standing, the body pushes through the feet and against the ground exactly as much as the ground pushes against the body through the feet. In dynamics, the balance is time-dependent and the solution to the system is through the equations of motion, which relate forces to motions. When the forces that produce the motions are fully known, the problem is characterized as forward dynamics, while when motions are known it is described as inverse dynamics.
In human bodies, the measurement of motions (via marker tracking or other means) is much less invasive and more feasible than measuring the muscle forces that produce the motions, so currently most analyses that depend on the musculoskeletal forces generated in human gait are typically treated as inverse dynamic problems (e.g. [72–74]). Further complicating the solution in multi-body problems (i.e. those with multiple, linked rigid bodies) is that the system can be interrogated from an external or internal perspective. The external view includes the action of the GRF to produce the movement of the body along a path, while the internal one includes both the joint reactions and the action of muscle forces in creating changes in joint angles. In mechanics, the relationship between the number of unknown variables and the number of known variables determines the difficulty of the solution. Fortunately, the internal and external perspectives can be exploited to develop solutions. For instance, in the single-stance foot, the GRF is the external force and produces the moment acting on the foot segment and reacted to by the ankle joint reaction force and moment (which are an internal force and moment). Once the ankle joint reaction force and moment are known, they can be combined with the motion of the shank to solve for the knee joint reaction and moment. The inverse dynamic solution allows for this time-dependent chain solution to determine the joint reactions.
2.3.4. Muscle redundancy solution
A system where the number of variables with unknown value equals the number of independent equations that fully define the system is determinant and a unique solution exists. For systems where there are more variables with unknown values than equations—called indeterminate—assumptions are made to provide more equations [75,76]. Many biological problems, including developing muscle forces, are indeterminate because there are more muscles capable of acting across a joint than necessary from a purely mechanical perspective, i.e. some muscles are (mechanically) redundant. If the desired result from an analysis that uses inverse dynamics is the development of muscle forces, then the joint moments determined from the inverse dynamic solution need to be apportioned to individual muscle forces. This is not a trivial undertaking, because the degree of redundancy is often quite large, so no general consensus on the most accurate criterion on which to base an apportionment algorithm currently exists [77].
Musculoskeletal modelling solutions typically include an algorithm that uses muscle parameters to apportion the joint moments to particular muscles (as muscle forces). Early solutions used PCSA to solve the muscle redundancy problem [75]: the apportionment was made based on the ratio of the PCSA of a particular muscle relative to the total PCSA of all muscles capable of acting (e.g. [61]). While effective from an algorithmic perspective, PCSA-based apportionment does not reflect the natural solution, because cross-sectional area approaches assume equal activation of all possible muscles.
Another algorithm for solving the muscle redundancy problem attempts to minimize the sum of muscle activations raised to a power [76]. The activation represents a proportion of total maximum strength the muscle can generate given the current state of the model (i.e. activation is generally bound between 0 and 1). If the muscle activation exponent is 1, i.e. a linear combination, force is allocated to the muscle with the greatest moment arm until it reaches its maximum allowed value, and then other muscles are recruited. Exponents greater than 1 will distribute forces more evenly across muscles available to produce the moment at a particular joint. The greater the power, the more even will be the distribution (PCSA-based activations are essentially an exponent of infinity). Other criteria for apportionment of the joint moments have also been advanced, such as electromyography (e.g. [78,79]), but activation-based methods remain the standard [76]. Until in vivo muscle forces can be reliably measured, all solutions to the muscle redundancy problem should be treated as hypotheses about how joints and muscles produce motion.
2.4. The output
A valuable aspect of musculoskeletal modelling (or of any model) is that variables of interest can be interrogated anywhere within the model. Standard output from a musculoskeletal model includes linear and angular displacements, velocities and accelerations for all body segments. Similarly, joint positions, velocities and accelerations are also tracked and can be exported for analysis. Muscle length, activation, moment arm and force are all possible outputs, as are joint reaction forces. Additionally, users can track the position of muscle origins/via points/insertions throughout the modelled motion. Conveniently, joint reaction forces can be output in the global coordinate system or the coordinate system of the relevant body segments.
Human locomotion requires the repetition of cycles of limb and body motions that can be reduced to a stride. For walking, a step is defined from an event (e.g. initial contact) on one foot to that same event on the contralateral foot and a stride is two contiguous steps. A stride can vary in duration, even at a constant velocity, so creating an average stride or evaluating the variability of stride parameters requires normalization of the strides to some consistent duration. The typical choice is per cent of stride cycle, where the stride duration is described in terms of equal portions of the entire stride (e.g. [80,81]). Evaluating the variability of parameters of interest (e.g. muscle force, joint reactions) can be an important step in validating the calibration of a solution. While some differences among trials of the same individual are expected, an evaluation of the variability in light of the question that the model is designed to answer should occur.
2.5. Human variability
A crucial advantage of musculoskeletal modelling is the opportunity to include aspects of morphological variation, including the production of subject-specific models. Model scaling and motion tracking incorporate some aspects of human variability, including differences in body mass, limb lengths/dimensions and kinematics, but other aspects are less obviously accessible for modification in musculoskeletal modelling software. These less accessible aspects, such as variation in bone or muscle morphology, are also important to consider. Quantifying human morphological (musculoskeletal) variation is at the heart of biological anthropology. For example, human (and hominin) pelvic morphology has a profound influence on bipedal capabilities [14,80], and has been the focus of intense scrutiny [41,81,82]. Humans are known to be sexually dimorphic in pelvic size and shape [83–85]. Furthermore, human pelvic morphology is also known to vary by population [86–88]. Sources of variation can be incorporated into musculoskeletal modelling to understand pathology [31] and differences in joint reaction forces [35]. For example, marker-informed parameter optimization and secondary scaling (described above) can change the width of the pelvis if experimental and virtual markers are included on bony landmarks of the pelvis, such as the ASISs. If the subject is wider than the generic model, the parameterization can widen the pelvis so that the virtual and experimental markers are more closely aligned. Scaling can also be accomplished manually either by applying factors to the three segment coordinate directions or by redefining the locations of specific anatomical landmarks associated with the segment. The former method moves all points proportionately while the latter allows for localized warping.
In addition to skeletal variation, aspects of muscle morphology that vary among humans may also be critical for accurate musculoskeletal modelling. For instance, variability in muscle parameters has been of concern for some time [89]. While repository models represent a particular person's morphology, people differ in important musculoskeletal modelling parameters such as muscle volume, optimal fibre length and pennation angle (table 1).
Table 1. Comparison of muscle parameters between the standard MoCap model in AMMR™ v. 2.3 (AnyBody Technology, Aalborg, Denmark; from Klein Horsman et al. [90] and Carbone et al. [34]) and Charles et al. [91]. Muscle names are per the MoCap AMMR v. 2.3. Where multiple muscle regions are present in the MoCap model, the values for the muscle volumes are summed while the optimal fibre lengths and pennation angles are averaged.Collapsemuscleoptimal fibre length (cm)pennation angle (°)muscle volume (ml)MoCapCharles—averageMoCapCharles—averageMoCapCharles—average
adductor brevis | 10.38 | 7.6 | 0.0 | 11 | 108.90 | 93 |
adductor longus | 10.57 | 11.0 | 0.0 | 12 | 159.56 | 159 |
adductor magnus | 9.97 | 23.1 | 0.0 | 12 | 559.25 | 567 |
biceps femoris caput breve | 9.14 | 10.9 | 0.0 | 9 | 107.95 | 92 |
biceps femoris caput longum | 8.54 | 20.4 | 29.9 | 11 | 232.01 | 194 |
extensor digitorum longus | 5.99 | 13.8 | 8.3 | 7 | 32.29 | 76 |
extensor hallucis longus | 5.98 | 10.6 | 14.4 | 7 | 36.27 | 21 |
flexor digitorum longus | 3.84 | 28.4 | 25.28 | |||
flexor hallucis longus | 5.20 | 30.1 | 161.50 | |||
gastrocnemius lateralis | 5.69 | 12.2 | 25.4 | 9 | 136.36 | 128 |
gastrocnemius medialis | 6.01 | 9.7 | 10.8 | 10 | 263.26 | 230 |
gemellus | 3.43 | 0.0 | 28.40 | |||
gluteus maximus | 13.57 | 0.0 | 936.55 | |||
gluteus medius | 6.76 | 5.3 | 674.71 | |||
gluteus minimus | 3.28 | 0.0 | 82.59 | |||
gracilis | 18.11 | 17.3 | 0.0 | 6 | 87.97 | 91 |
iliacus | 8.16 | 8.8 | 203.13 | |||
obturator externus inferior | 6.88 | 0.0 | 37.88 | |||
obturator externus superior | 2.77 | 0.0 | 68.18 | |||
obturator internus | 2.05 | 0.0 | 52.08 | |||
pectineus | 9.50 | 0.0 | 64.58 | |||
peroneus brevis | 4.54 | 23.1 | 86.09 | |||
peroneus longus | 5.08 | 15.8 | 121.52 | |||
peroneus tertius | 4.27 | 19.1 | 26.52 | |||
piriformis | 3.88 | 0.0 | 31.25 | |||
plantaris | 4.77 | 0.0 | 11.36 | |||
popliteus | 2.40 | 7.4 | 0.0 | 8 | 25.57 | 15 |
psoas major | 9.92 | 13.4 | 193.18 | |||
psoas minor | 5.91 | 0.0 | 6.63 | |||
quadratus femoris | 3.37 | 0.0 | 49.24 | |||
rectus femoris | 7.84 | 14.2 | 21.9 | 8 | 226.33 | 249 |
sartorius | 34.71 | 40.8 | 0.0 | 205.49 | 143 | |
semimembranosus | 8.09 | 15.8 | 25.0 | 12 | 138.26 | 244 |
semitendinosus | 14.16 | 18.3 | 0.0 | 8 | 208.33 | 186 |
soleus | 4.40 | 14.6 | 61.6 | 12 | 792.91 | 461 |
tensor fasciae latae | 9.49 | 0.0 | 83.33 | |||
tibialis anterior | 6.83 | 13.7 | 9.6 | 7 | 181.49 | 129 |
tibialis posterior | 3.78 | 34.2 | 163.36 | |||
vastus intermedius | 7.68 | 18.1 | 11.8 | 11 | 292.61 | 521 |
vastus lateralis | 6.69 | 19.6 | 0.0 | 15 | 583.33 | 606 |
vastus medialis | 7.16 | 15.9 | 0.0 | 14 | 454.27 | 415 |
Several possibilities for future work that explores the impact of morphological variation on muscle and joint forces are worth noting. For instance, while some anthropological work has elucidated the interactions of morphology and walking conditions (e.g. velocity, changes of direction, burdens and gradients) on energy expenditure (e.g. [92–94]), much less is known about how morphology and walking condition factors interact to affect joint and muscle forces (although see [31] for an example). People vary morphologically and walking under varied conditions is an activity of daily living, so this area is an important one for many disciplines to pursue including biological anthropology, biomedical engineering and product design. Perhaps more importantly, collecting the requisite data to inform a musculoskeletal model is time consuming, invasive and dependent on (at this time) fragile laboratory equipment. Consequently, acquiring samples that represent the breadth of human morphological variation is unlikely to occur in the near future. Once validated, though, a musculoskeletal model can be modified to assess the impact of these known morphological differences on muscle and joint forces. Extending that line of research further, a longer-term goal would be use musculoskeletal modelling to understand the influence of the morphological differences between modern humans and extinct hominins on muscle and joint forces, an area that has sought biomechanical solutions without effective software support for decades [12,14,61,95,96].
2.6. Current study
Here, we demonstrate the process and output of one human subject during walking using AnyBody Technology's musculoskeletal modelling software (Anybody Modeling System™; AnyBody Technology, Aalborg, Denmark). We present several of the possible outputs that could be of interest to biological anthropologists. This includes joint kinematics, joint reaction forces and a selection of lower limb muscle force magnitudes. We make some recommendations for future researchers to help facilitate their use of software and output. Additionally, by presenting a set of possible outputs, we demonstrate the motivation for undertaking musculoskeletal modelling and that results may have a direct bearing on anthropological questions, or as intermediate steps to other analyses (e.g. finite-element analysis). Additionally, while engineering has much to offer biological anthropologists, the latter have an appreciation for human variation which may enhance future generations of musculoskeletal models.
3. Material and methods
3.1. Participant
One healthy participant was used for this demonstration, but he is part of a larger study aimed at examining human walking. The participant was 36 years old (male; body mass, 74.3 kg; stature, 172 cm) and reported being free from lower-limb injuries. The University of Washington's Institutional Review Board approved all aspects of this study (IRB no. STUDY00001125).
3.2. Experimental protocol
Kinetic and kinematic data were measured using a four-force plate (Kistler, Switzerland), 10-camera motion capture system (Qualisys, Sweden) in the Amplifying Movement & Performance (AMP) laboratory of the University of Washington. Thirty infrared-retroreflective 5 mm hemispherical markers were affixed to anatomical landmarks and used to track motion through the laboratory (figure 1; electronic supplementary material). The participant walked unshod the length (10 m) of the gait laboratory five times at his self-selected preferred walking speed (1.15 m s−1). The participant walked in a straight path with his eyes directed at a point on the far wall of the laboratory and contacted the surface of each force plate with a single foot. Trials with multiple feet on a single force plate or ones where the stance foot crossed plates were immediately redone. The participant was allowed several attempts prior to data collection to become familiar with the protocol. All trials exhibited similar kinematics and kinetics (data not shown), so we selected one trial to use in this demonstration.
Marker data and GRFs were collected at 120 Hz and 1200 Hz, respectively, and then were filtered at 5 Hz and 10 Hz, respectively, with a fourth-order, low-pass zero-lag Butterworth filter. Calibration of the system yielded a limitation in its fidelity for marker data of 1 mm and force data of ±2.5 N for the direction of travel (X), ±5 N side (Y) and ±25 N vertical (Z). All data were exported from the Qualisys Track Manager software in C3D file format, which can be read directly by the modelling software (provided in electronic supplementary material).
3.3. Musculoskeletal model
We used the ‘MoCap’ model from the AnyBody Managed Model Repository™ (AMMR v. 2.3; AnyBody Technology, Aalborg, Denmark; http://www.anybodytech.com/software/model-repository-ammr), which is a multi-trial, full-body, motion-capture-driven human gait model, and the commercially available AnyBody Modelling System™ (v. 7.3; AnyBody Technology, Aalborg, Denmark) to calculate forces in the lower limb [34,97]. The MoCap model that we used for this demonstration has been used to assess human gait in numerous applications (e.g. [98,99]). The MoCap model is assembled from a trunk and left and right lower limb components. The trunk includes a lumbar region, trunk, neck and head. The lower limbs consist of the following segments: thigh, patella, shank, talus and foot. Each lower limb has six total DOFs, including all three rotations at the hip and one each at the knee (flexion/extension), ankle (plantar flexion/dorsiflexion) and subtalar (inversion/eversion) joint. The pelvis (relative to the ground) has six DOFs, three translational and three rotational. Forty-one anatomically defined, bilateral muscles (table 1) are represented by 169 muscle elements in each model lower limb (e.g. gluteus medius is composed of 12 separate muscle element actuators) [34]. The muscle elements in this model generate force as a function of their activation and a strength parameter which is appropriate for human walking [97,100]. We defined a virtual motion-tracking marker for each of the experimentally measured motion markers (electronic supplementary material). The polynomial algorithm that was employed to solve the muscle redundancy problem apportioned force to muscles using an objective function with power = 3. (For more information regarding the model, see §5.3.)
3.4. Musculoskeletal simulation
The first step of the simulation was to scale the model to match the dimensions of the subject while simultaneously optimizing the marker coordinate and joint rotation axes' locations [68]. After initial scaling, the values of these parameters (i.e. segment dimensions, marker locations, rotation axes' locations) were determined using the marker motion data from a trial. Values for all parameters were optimized to minimize the sum of square distances between model marker positions and experimental marker positions. Following convergence of parameter optimization, we conducted an inverse kinematics analysis. Finally, we conducted the inverse dynamic analysis, which includes apportioning muscle forces using the algorithm that minimizes the cost function as a sum of muscle activation values raised to the third power. This value has been shown to allocate muscle activations in a way that reflects experimental data [77,100,101].
3.5. Output variables and data processing
We output several kinematic and kinetic variables, including joint kinematics, filtered GRF components and joint reaction forces (electronic supplementary material). Additionally, we output the force magnitude for a selection of the muscle elements for the left lower limb during an entire stride cycle of walking. Muscle element forces within anatomically defined muscles were summed for presentation (e.g. 12 elements that represent gluteus medius). All time-related curves were resampled to 1% increments of the stride cycle.
4. Results
GRF component curves across the gait cycle are presented in figure 2. As is typical for human walking, the vertical component has two peaks (early stance peak, 719 N; late stance peak, 778 N) on either side of a valley (633 N) which occurs around the middle of stance phase. The average of the two peaks is 103% of body weight (729 N), while the valley represents 87% of body weight.
4. 결과
보행 주기에 걸친 GRF 구성 요소 곡선은 그림 2 에 나와 있습니다. 사람의 보행에서 흔히 볼 수 있는 것처럼, 수직 성분은 스탠스 단계의 중간쯤에 발생하는 계곡(633 N)의 양쪽에 두 개의 피크(초기 스탠스 피크, 719 N, 후기 스탠스 피크, 778 N)가 있습니다. 두 피크의 평균은 체중의 103%(729N)이며, 계곡은 체중의 87%를 나타냅니다.
Figure 2. Filtered ground reaction force components. Grey vertical line represents moment of toe-off.
Pelvic rotation and obliquity and hip, knee and ankle kinematics are provided in figure 3. All kinematics are of the left lower limb throughout a single stride. Pelvic motion is described in terms of the right hip (previous stance limb) relative to the left hip (current stance limb) [102]. Hip, knee and ankle joint reaction forces, which include contributions from muscle forces, are provided in figure 4. In each joint, the vertical joint reaction force is greater in magnitude than the anteroposterior and mediolateral components. All three joints show a large peak during the second half of the stance phase, and this late peak is significantly larger in the ankle. Muscle force profiles for 18 anatomically defined muscles across the stride cycle are provided in figure 5. During the stance phase of walking, gastrocnemius, soleus and gluteus medius generate the largest forces. All output data are provided in the electronic supplementary material.
Figure 3. Angular kinematics for the pelvic segment and hip, knee and ankle joints. Vertical grey line represents toe-off.
Figure 4. Joint reaction forces expressed in the global coordinate system. Vertical grey line represents toe-off.
Figure 5. Muscle forces for 18 anatomically defined muscles.
5. Discussion
5.1. Model output
The subject that we used to demonstrate the musculoskeletal modelling approach to modelling walking is an adult male without gait pathology and within species-typical anthropometric characteristics. Consequently, the model-derived motions should accord with those determined experimentally and extensively described in the last 40 years of gait research (e.g. [66,103–105]). The model results presented above do meet this requirement, as discussed below in detail. We focus on the model-derived kinematics, because these kinematics are the foundation for all subsequent analyses in musculoskeletal models. Any analysis should begin by verifying that model-derived variables align with those obtained empirically, and the cause of any discrepancies should be understood or noted as a potential limitation. We belabour the details in the following paragraphs in order to demonstrate the basic approach to model validation.
In our subject at heel strike, the pelvis is rotated approximately 10° in the horizontal plane such that the right hip is posterior to the left hip (e.g. the current stance limb). The pelvis then swings through a neutral position and, while leading up to the next right heel strike, into an orientation with the right hip anterior to left hip (approx. 10°). The pelvis drops in the frontal plane with the right hip inferior to the left (stance) hip soon after heel strike as the right lower limb moves into the swing phase with pelvic obliquity ranging between ± 4°. The pelvis remains in this position until the subsequent right heel strike, after which the pelvis approaches the neutral position. Model-derived pelvic rotation and obliquity are similar to empirically determined values (e.g. [103]).
In the sagittal plane, the left hip begins the stance phase in flexion of approximately 17° and moves through a neutral position and into extension of approximately 17°. Winter [66] describes a less symmetrical flexion–extension cycle at natural cadence (mean values of approx. 20° of flexion and 11° of extension), but the model-derived hip rotations are well within one standard deviation of his mean values. During most of the stance phase, the left thigh is adducted less than 6° in accord with others (e.g. [103]).
While the pelvic and hip joints have three rotational DOFs, the knee and the ankle only rotate in the sagittal plane. The left knee flexes slightly (10°) early in the stance phase, but maintains a near fully extended position through midstance, after which it begins to flex in preparation for the swing phase where flexion exceeds 60°. After the initial heel strike and as the forefoot contacts the substrate, the ankle plantar flexes slightly (less than 5°) and then assumes a dorsiflexed position for most of stance. Rapid plantar flexion occurs in late stance as the body is propelled forwards, reaching a plantar flexion value of greater than 10°. The ankle returns to a dorsiflexed position for swing. These values for both knee and ankle angles are typical (e.g. [66,103,105]).
Consequently, the kinematics of the model match those that we expect empirically. While this check is not sufficient to ensure that the internal (muscle and joint reaction forces) are reasonable, it is a necessary condition and validation for any model. Incorrect kinematics guarantee incorrect or suspect muscle forces.
Joint reaction forces have been empirically measured only with instrumented joint implants. By definition, the surgery creates conditions that limit applicability of the results, but no other method currently exists to measure force. Pedersen et al. [45] measured hip joint contact forces in a 72-year-old man (63.2 kg) with a hip replacement and found that the vertical force was 1750 N, considerably less than the almost 3500 N we found with the musculoskeletal model. Important differences between this patient and our subject exist, including age (72 versus 36 years), body mass (62.3 versus 74.3 kg), health status (58 days post hip replacement surgery versus healthy) and walking speed (0.89 versus 1.15 m s−1). While hip and knee joint reaction forces presented here are larger than in vivo values using instrumented prostheses [77,106,107], our model-derived joint forces are, however, similar to other models of healthy adults [35,107].
Although it is difficult to compare muscle forces with experimental data, it is possible to consider the pattern of activation with reference to the muscle's action. The muscles that generate the highest force in the model are gluteus medius, soleus and gastrocnemius. Gluteus medius plays a critical role in bipedal walking by preventing pelvic drop during single stance [103,105]. The gluteus medius muscle profile from the model is double peaked, and the timing of the two peaks is coincident with the two peaks in the vertical GRF, indicating the muscle's support of the pelvis during single stance. Gastrocnemius and soleus plantar flex the ankle; when the foot is in contact with the ground and the hip is extended, plantar flexing the ankle moves the body anteriorly [103,105]. The gastrocnemius and soleus muscles show the greatest force generation during the propulsive portion (second half) of the stance phase, in keeping with our expectation from their function.
Other muscles also display moderately high forces. Tibialis anterior, which dorsiflexes the ankle, shows muscle force generation immediately after heel strike when the body begins to move over the foot, and then is mostly inactive until the second half of the swing phase. The quadriceps demonstrate a force peak early in stance phase, and again at the stance–swing transition, potentially relating to their role in extending and stabilizing the knee [105]. The hamstring muscles show the greatest peak at the end of swing phase and into early stance phase. Activation of these muscles at the end of swing phase is associated with decelerating the limb prior to heel strike. Early in stance the GRF is directed anterior to the hip and knee, and the hamstring muscles are active to prevent further flexion (hip) and extension (knee) of these joints. The hip flexors are active late in the stance phase, and may be active to counteract hamstring hip extension moment as the hamstrings flex the knee in preparation for swing phase. The adductors show initial peak in stance (adductor magnus) and late stance (adductor longus), but generally smaller than forces in other muscles.
5.2. Input data choices
The model-derived motions and joint reaction and muscle forces align with our expectations from the experimental data in the literature. We achieved these results through careful attention to detail when we collected the data we used as inputs to the musculoskeletal model. We provide below some brief suggestions from our experience. Because motions are derived from the change in position of these markers, marker placement on the segment should be chosen to define fully the motion of the segment and to reduce marker positional noise as much as possible. Positional noise arises from the inherent error of measurement of the marker's position in space (system error) and from movement of the surface (skin or clothing), where the markers are placed relative to the underlying rigid segment (e.g. [108–110]). System error is reduced through calibration; for our system, it is less than 1 mm for any marker. Movement between the external (skin) surface and the bone is more complicated to reduce, but may not be clinically relevant [110].
Marker position is a compromise with soft tissue coverage of underlying bone, landmark palpability (to aid placement) and distance between markers (such that marker positional error is much less than the distance between markers) all being important. The shank and foot have palpable landmarks that are covered by minimal soft tissue (e.g. the malleoli), but the thigh (e.g. greater trochanter) and, especially, the pelvis are more challenging segments on which to affix markers. This problem is exacerbated when the subjects of interest have more soft tissue (muscle or fat) overlying the bony landmark, and considerable soft tissue coverage is common in the hip–pelvis region. One approach to ameliorate some of this difficulty is to use a marker plate which has four (or more) markers mounted on a hard plastic plate that is then strapped to the segment. If a plate can be affixed firmly to the segment, as it can be to the thigh, the marker plate eliminates between-marker error. Unfortunately, marker plates are difficult to affix firmly to the pelvis. Marker placement issues, consequently, are most difficult to overcome to track the motion of the pelvis. Yet, for musculoskeletal models of human locomotion, pelvic motions are critical components of system performance. We achieved the motions shown in figure 3 through extensive protocol development and testing.
Another critical input to the simulation are the GRFs. It is possible to conduct an inverse dynamics solution of the single support phase of walking with a single force plate (e.g. [111]), but extending the analysis to the double support portions of the stride is more difficult with one force plate because it would require estimation of the unknown contribution of the other stance foot. Consequently, any research that requires information about double support should optimally be conducted with multiple force plates. The placement of the force plates relative to each other should be appropriate to the population and gait. Stride lengths for adults or running gaits are longer than for children or walking, and plates should be positioned to prevent subjects ‘targeting’ foot placement or large numbers of repeated trials. Data analysis is less complicated if stance occurs on one plate, but the data from two immediately adjacent plates can be used as long as there is no gap between them [112]. Multiple feet contacting a single plate is also possible to resolve with additional information (such as localized pressure sensors) or predictive methods (e.g. [113,114]). All of these remediation techniques are time-intensive. The laboratory in which the presented data were collected has four floor-embedded force plates that are positioned to allow normal stride lengths for adults when walking. The subject and trial that we used for this demonstration is part of a much larger study [115]. Nonetheless, we assessed each foot placement on all four plates for each trial before moving on to the next trial. Such data checking is a tedious effort at the time of data collection, but well worth the effort to ensure all data are suitable for future analyses. Using four force plates, we had data frames before and after the data that we used in the analysis. This allowed us to discard the initiation and termination frames of the simulation, which are frequently unusable. Consequently, researchers are advised to locate force plates well and test the placement prior to data collection.
Finally, researchers may consider collecting complete anthropometrics of subjects and using these in the initial scaling as the starting points for the parameterization. People vary considerably in terms of their proportions [35,37], so standard models, such as the model we used here, represent a morphology that may not match the dimensions of any particular subject. If your model starts with segment lengths that are substantially different from the segment lengths of the model, the parameterization algorithm may not be able to converge on a suitable solution. Two failures can occur: the parameterization can fail to converge, which produces an obvious error message, or the parameterization converges, but, in order to converge, the simulation requires changes to the morphology or motions that are inconsistent with the experimental data or subject dimensions. This latter failure can be insidious. An absolutely critical step in model validation is to check the motions and shapes of the segments after parameterization. A poorly parameterized model is incorrect and will not produce acceptable motion or joint and muscle forces.
5.3. Modelling choices
In order to produce this demonstration, we made a number of choices with regard to the model which bear discussion. The most important of these was the choice to use a standard, repository model, in this case the MoCap model of AMMR v. 2.3. By using the MoCap model we accepted many assumptions and compromises, but gained the experience and efforts of many dedicated musculoskeletal modellers. Nonetheless, the responsibility of assessing the ability of any model to answer the question of interest rests with the researchers using the model. Our intention was to demonstrate some of the decisions and validations needed in musculoskeletal modelling of human walking and our choices reflect that.
The first of these choices is the theory on which the muscle actuation rests. The muscle models used here generate force as a function of maximum force and activation. More sophisticated muscle models, such as the Hill muscle model, are available. The Hill muscle models are thought to provide the most complete model of muscle behaviour [116], but the Hill model requires several parameters for its full definition, including: maximum isometric force, tendon slack length, optimal fibre length, pennation angle at optimal fibre length, activation and deactivation times, as well as maximum contraction velocity. The number of parameters that must be defined for a Hill muscle model, as well as the sensitivity of muscle behaviour to parameter values, means that considerable effort must be expended to parameterize and optimize Hill-type models. As with all models, poorly informed inputs (i.e. muscle parameters) lead to unstable results. This means that for certain activities it may be advisable to use a simpler, robust model rather than a more complex, unreliable one. For human activities during which muscles operate at slow contraction velocities and near their optimal fibre length (e.g. human walking), simple muscle models that depend only on geometry and a maximum force parameter are preferable [100]. Consequently, we accepted the use of the simple muscle model employed in the MoCap gait model.
Another choice was the selection of the scaling method. Understanding which scaling method to use (e.g. mass–length–fat), requires understanding the underlying assumptions of the method (e.g. distribution of mass to segments), to which model parameters the scaling applies (e.g. which segment length scale to which muscle length) and how the model uses the scaling to define parameters (e.g. are all parameters scaled or just a subset?). We made many other choices by not changing components that are inherent to the generic MoCap model, including acceptance of the form of the algorithm used to solve the muscle redundancy problem. Such choices are often embedded in standardized models, and the user does not have to modify any of them in order for the model to run to completion. This makes ensuring that the question of interest can be addressed by the model become an intentional assessment. For instance, the standard MoCap model has only one rotational DOF for the knee, which works well to assess overall walking parameters but is less well suited to a detailed study of the knee. Furthermore, standardized models are improved as more experience and data become available (e.g. [97]). For instance, the MoCap model's leg (TLEM v. 2.1) has been modified to reflect muscle paths more accurately by the inclusion of wrapping surfaces.
We also relied on the standardized definition of available outcome variables such as muscle forces and joint motions. Given that most software allows for forces and motions to be expressed in different ways (at a minimum, in different coordinate systems) it is imperative to understand the context of the reported variable. Is an output reported in the global or segment coordinate system? Is a joint force an internal force or a reaction? In what direction does the muscle force act? As with all big data, extensive data exploration and a complete understanding of the theories (mechanics), procedures (inverse dynamics) and systems (human anatomy) is necessary to draw conclusions from a musculoskeletal modelling simulation.
5.4. The role of anthropology in musculoskeletal modelling
Finally, it is important to make clear that in this demonstration we have used a walking trial of one subject and focused on the subject's left side. The kinematics and GRFs produced by this individual's right side and by other walking trials demonstrate a consistent pattern, but were not numerically identical to those reproduced here (electronic supplementary material). Consequently, the muscle forces and joint reactions are also somewhat different. The degree of intraindividual variability in gait parameters is influenced by such characteristics as age and health, but, even for this healthy adult male, peak forces within a muscle can vary. In muscles that exert relatively small forces (less than 500 N), step-to-step values can vary by more than a factor of 2 (electronic supplementary material). Muscles which show greater force generations (e.g. gastrocnemius; greater than 1000 N), step-to-step variability is lower (a factor of approx. 1–1.4; electronic supplementary material). This serves to underscore the importance of collecting and analysing multiple steps when the intent is to report patterns and forces (not the intent here).
While scaling addresses some aspects of variability and improves model predictions [35,69], scaling itself does not solve all the problems. Interindividual variability in such important variables as muscle parameters has been appreciated for over a decade [89], but the consequences of such variability for musculoskeletal modelling remains unclear. For instance, muscle strength as implemented in the solution that we use for this report depends on three factors: muscle volume, optimal fibre length and pennation angle. The MoCap model, which we used in our demonstration, uses a set of values derived from several sources but is critically informed by Klein Horsman et al. [90] and Carbone et al. [34]. These data are from measurements of cadavers who were 77 and 85 years old at the time of death, respectively. Cadaveric data are widely used in musculoskeletal models (e.g. [53,99,117]), but muscle parameters from cadavers may not represent those from individuals who are not at imminent risk of death. Table 1 includes values from Klein Horsman et al. [34] and Carbone et al. [34] and from a recent analysis of 10 healthy individuals who were evaluated with diffusion tensor imaging [91]. Gastrocnemius and soleus, the principal drivers of forward propulsion, vary in all three factors. Importantly, the complexity of gait mechanics that makes musculoskeletal modelling an attractive solution also makes understanding the effect of this variation impossible without simulating the difference. Larger volume, smaller pennation angle and shorter optimal fibre lengths all act to increase muscle strength, thereby increasing the allocation of force to it by the algorithm. But, the algorithm uses a muscle's strength relative to all other muscles' strengths, making it impossible to predict the impact of muscle parameter variation on muscle forces or joint reactions. Variation in muscle parameters is but one of the many ways that people vary in their musculoskeletal anatomy. Leveraging the expertise of biological anthropologists to assess the impact of human biological diversity on musculoskeletal modelling can only be beneficial for the entire community of researchers.
6. Conclusion
Although motivated by different specific goals, biological anthropologists and engineers (principally biomedical) both pursue questions that revolve around understanding the internal forces of the primate musculoskeletal system. Engineering questions tend to be species-specific, focusing on living humans. Biological anthropologists are also interested in living humans, but frequently ask questions that include greater taxonomic diversity. Research in both fields can be specific to the details of the individual (e.g. patient-specific orthopaedics, particular hominin fossils) or generic representatives of groups (e.g. sex-specific running shoes, comparisons of primate taxa).
Musculoskeletal modelling is powerful because of the flexibility it provides at every stage of research. Building a new model requires intensive effort, but may be necessary based on the question. Previously validated models from repositories lower barriers to modelling experiments, but researchers should be aware of choices made by model creators. Like all modelling exercises, serious and intentional thought must be dedicated to the multitude of modelling choices. This is not to say that all models need to be complex. If a simpler model adequately describes the phenomenon of interest and answers the question at hand, then a simpler model may be preferable.
Critical to support rigorous simulations, musculoskeletal models include changeable parameters that describe all aspects of the musculoskeletal system. This is a purposeful feature included by engineers that drove musculoskeletal modelling and an attractive feature for biological anthropologists. Model users can change the size and shape of the skeletal elements, dimensions and inertial properties of body segments, DOFs and range of motion of joints, action of passive soft tissue structures, the geometry and force capacity of muscles, as well as numerous other parameters. Generic models from repositories (such as the one used here) can be scaled and reshaped to capture many aspects of human variation (e.g. limb dimensions, mass), and hence are well suited to many questions. More effort can be dedicated to refine aspects of a model that are of particular interest (e.g. change the morphology of the bony pelvis). Modelling flexibility naturally extends to the information that can be derived from simulations. Forces generated by individual muscles (or portions of muscles) as well as internal to joints are standard output, as is the ability to express these quantities in different coordinate systems. Until internal forces of the body can be measured non-invasively, musculoskeletal modelling will remain the pre-eminent approach to these questions.
Future collaborations between engineers and biological anthropologists to create better models of human movement have the potential to expand the understanding of mobility in both fields. Subject-specific models are internally consistent but not generally representative of people, so extrapolating should only be done with this limitation in mind. Developing a model from scratch for a fundamentally different organism (e.g. a chimpanzee) is a significant investment.
Data accessibility
Additional information is available as the electronic supplementary material. Raw data and simulation output are available at the assigned figshare link.
Authors' contributions
A.D.S. and P.A.K. conceived the project. S.G.L. and P.A.K. designed and carried out the data collection. A.D.S. and P.A.K. carried out the simulations and data analysis. A.D.S. and P.A.K. wrote the draft and all authors contributed to the final version of the manuscript.
Competing interests
We declare we have no competing interests.
Funding
P.A.K. and S.G.L. thank S.K. Benirschke, MD, Jerome Debs Endowment Chair in Orthopaedic Traumatology at the University of Washington, for the ongoing support of this research and for our many conversations about foot form and function.
Acknowledgements
The authors would like to thank the participants of this study for their time and patience during the protocol. The authors would also like to thank the reviewers and editors for comments that improved this manuscript.
Footnotes
One contribution of 12 to a theme issue ‘Biological anthroengineering’.
Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5486386.
Review articles
A review of musculoskeletal modelling of human locomotion
,
and
Published:13 August 2021https://doi.org/10.1098/rsfs.2020.0060
Abstract
Locomotion through the environment is important because movement provides access to key resources, including food, shelter and mates. Central to many locomotion-focused questions is the need to understand internal forces, particularly muscle forces and joint reactions. Musculoskeletal modelling, which typically harnesses the power of inverse dynamics, unites experimental data that are collected on living subjects with virtual models of their morphology. The inputs required for producing good musculoskeletal models include body geometry, muscle parameters, motion variables and ground reaction forces. This methodological approach is critically informed by both biological anthropology, with its focus on variation in human form and function, and mechanical engineering, with a focus on the application of Newtonian mechanics to current problems. Here, we demonstrate the application of a musculoskeletal modelling approach to human walking using the data of a single male subject. Furthermore, we discuss the decisions required to build the model, including how to customize the musculoskeletal model, and suggest cautions that both biological anthropologists and engineers who are interested in this topic should consider.
1. Introduction
A principal research focus within biological anthropology is interpreting skeletal variation within the context of behavioural diversity, including variation in diet, disease processes and activity patterns among living and extinct primates, including humans and hominins [1–4]. Locomotor behaviour and musculoskeletal morphology have been a central focus of this area of inquiry because locomotion provides primates with access to key evolutionary resources including food, shelter and mates. Furthermore, bipedalism is considered a defining character of the hominin lineage, cementing its importance within the field [5–7]. The challenge of elucidating form–behaviour relationships are the complexities of the intervening and underlying functions (e.g. joint range of motion) and biomechanics (e.g. joint reaction forces).
Early anthropological research connecting form and function was both typological and qualitative in nature; however, anthropologists have increasingly adopted mechanical engineering approaches to answer questions about primate locomotion generally, and hominin locomotion in particular (see [4] for a historical review). Within anthropology, two major subfields of mechanics have been applied to locomotor systems: Newtonian mechanics and solid continuum mechanics. Newtonian mechanics is concerned with the description of the motion of solid bodies under the influence of a system of forces, while solid continuum mechanics quantifies the behaviour (e.g. deformation) of solid materials when subjected to forces [8,9]. Lovejoy and colleagues [10–13] were early adopters within anthropology of both Newtonian mechanics and solid continuum mechanics [4]. Newtonian mechanics is routinely used by anthropologists to investigate a variety of questions (e.g. [14–19]). Beam theory, an application of solid continuum mechanics, has been used to estimate the structural capacity of long bones and is regularly exploited by anthropologists (see [4] for a review) to reconstruct hominin locomotor adaptations (e.g. [13–17]) as well as modern human mobility patterns (e.g. [18–24]).
The application of Newtonian and continuum mechanics is, of course, a core competency of mechanical and civil structural engineering, and the methods that anthropologists use derive directly from these engineering fields. The techniques and theories that underlie their application were developed, however, to design and analyse human-made products using inert materials. While such structures and machines tend to be faithfully produced based on design specifications, humans are not clones. Rather, each individual is completed from a standard template that has, as an inherent feature, the ability to respond to the environment. This means that whereas human-made products generally have low tolerance for variance, among humans, variation is the rule. Furthermore, the ‘standard’ template itself varies among ancestral lineages, populations and species because of the evolutionary history of each (e.g. [25–27]). The biomedical engineering world has responded to this with the use of patient-specific models because, in patients, a direct relationship exists between a particular morphology (i.e. that of the patient) and the question to be answered (e.g. which implant configuration limits strain on the bone?) with the musculoskeletal model (e.g. [28–35]). Not all questions, however, concern specific individuals. In anthropology, the focus may be comparing species (or populations) to understand responses to selective pressures. In engineering, running shoes, for example, are not designed for particular individuals, who vary in foot function [36], but rather for people generally. And people are not just scaled versions of each other [37]. Thus, an appreciation of human variation, a core competency of biological anthropology, is a critical component of any modelling exercise.
At the heart of many locomotor-focused questions, both engineering and anthropological, is a need to evaluate internal (e.g. muscle, joint contact) forces. Finite-element analysis of skeletal elements [38,39], estimating locomotor costs [40,41] and evaluating task performance [42,43] all require muscle forces to be known or estimated. Measuring muscle and joint contact forces in living animals is, however, exceptionally challenging because it requires surgical intervention to implant the relevant sensors into the body [44,45]. As an additional obstacle, these techniques are normally limited to a single joint or muscle/tendon within a single study (see [44] for examples). Consequently, internal forces must be estimated. Poorly estimated loads are likely to produce misleading results, interpretations and conclusions. Consequently, estimating internal forces (e.g. muscle forces) in humans is of critical interest to both anthropologists and engineers.
1.1. Musculoskeletal modelling: a path forwards
Originally driven by a clinical interest in gait pathologies, musculoskeletal modelling has emerged as the pre-eminent technique to elucidate the internal details of human (and other animal) movement [46]. These models build upon dynamic analyses of linked rigid-segment models ([47] and described below). Rigid-segment models represent the body as a series of solid body segments which are linked together at joints. Musculoskeletal modelling extends rigid-segment models by including detailed models of individual muscles, as well as models of neurological control. In extending the dynamic approach, musculoskeletal modelling solves the muscle redundancy problem based on muscle parameters under a specified minimization criterion [48–51]. The muscle redundancy problem refers to the fact that there are more muscles than necessary to actuate movements at joints [48,49]. In solving for muscle forces, this leads to an indeterminate system of equations in which there are more unknowns (muscle forces) than constraints (equations) [51], thus requiring optimization solutions. Results of musculoskeletal modelling can reveal individual sequences of muscle excitation and activation, as well as estimates of muscle forces and bone-to-bone contact forces.
In a series of papers, Pedersen et al. (see [45]) developed an early version of this type of model and used it to estimate lower limb muscle forces and hip joint contact forces. The hip joint contact forces were validated using an instrumented hip joint replacement in a single elderly subject [45]. Consequently, these data provided an early demonstration of the validity of musculoskeletal modelling as an approach to determining internal forces. Since then musculoskeletal modelling has been used for numerous projects exploring human motion and locomotion [31,43,52–56] as well as that of non-human primates [57–59] and early hominins [60]. Wang et al. [61] applied this technique to the hominin fossil record, but their approach was necessarily limited by the lack of software available to carry out the complex calculations. Fortunately, several software packages (e.g. AnyBody Modeling System™, OpenSim, FreeBody) are now available to carry out musculoskeletal modelling [45,62,63].
Accurate musculoskeletal models demand both understanding of the principles of Newtonian mechanics and an astute appreciation of anatomy and anatomical variation. Thus, successful musculoskeletal modelling is clearly situated at the interface between engineering (i.e. mechanics) and biological anthropology (i.e. human variation). Here, we provide an overview of this effort from an anthroengineering perspective, blending knowledge from both domains with a special emphasis on a task—understanding human walking—that is of equal interest to both groups.
2. Background
2.1. The model
In musculoskeletal modelling used to solve for muscle forces, the body of interest (e.g. the whole human body) is modelled as a series of rigid segments connected by joints and actuated by muscles [34,46] (figure 1). This body moves in a virtual space whose dimensions are defined by a global Cartesian coordinate system. Models can include as few or as many body segments as necessary to answer the question of interest. For walking, models minimally include segments for the pelvis (and trunk/arms/head, often aggregated) as well as left and right feet, shanks and thighs. More complete models often also include segments that represent the patella and talus, as well as a head–arm–trunk (HAT) segment separate from the pelvic segment. Additionally, the upper part of the body may also be modelled to include segments for the lumbar region, thorax, arm, forearm, hand, neck and head segments. Body segments are defined in terms of their physical dimensions, mass, location of centre of mass and moments of inertia. Each body segment also includes its own Cartesian coordinate system, usually aligned with standard anatomical directions of the segment (e.g. proximodistal, anteroposterior, mediolateral). Other coordinate systems can be defined relative to that of the segment or to the global system and allow for ease in defining model geometry. The segments further rely on the assumption that they do not deform during locomotion (i.e. are rigid). In the software modelling/simulation environment, body segments are almost universally represented by virtual surface models of the relevant skeletal element(s).
Figure 1. MoCap model mid-simulation of walking. The image shows the model (rigid body segments represented by skeletal elements and visual representation of muscle models), four force plates (grey), the ground reaction force vector (blue line through model), force plate readings (blue lines below the force plate), experimental marker locations (blue spheres), virtual markers (red spheres with coordinate system arrows) and global coordinate system (yellow arrows).
Body segments are connected to other segments at joints, which represent articulations (e.g. hip, knee, ankle) within the body. Joints are frictionless and non-deformable and, consequently, forces and moments transfer across these joints in an equal and opposite manner. A joint definition includes its location relative to the body segments it connects, as well as the degrees of freedom (DOFs) permissible at the joint (one, two or three axes of rotations and, less commonly, translations). The direction of rotation axes are usually expressed in terms of a local (body segment) coordinate system.
With the body segments and joints defined, the skeletal portion of the musculoskeletal model is sufficient to conduct dynamic analyses. Musculoskeletal models extend the analysis by including models of individual muscles or even portions of muscles (figure 1). Muscle models range in their complexity and, as with other parts of the model, the necessary complexity depends on the motion being simulated. (For a description of how muscle complexity can be modelled, see §5.3.) A muscle model minimally requires geometry (i.e. path definition) and a parameter representing its maximum possible force. More complex muscle models also include parameters for pennation angle, optimal fibre length, physiological cross-sectional area (PCSA) and activation and deactivation times. Muscle geometry must include the coordinates of its origin (location relative to the proximal segment) and insertion (location relative to the distal segment), but muscles can be geometrically more complex than a straight line path between origin and insertion. The most basic of these complexities derives from a muscle with a broad origin (e.g. gluteus medius) or insertion (e.g. adductor magnus). In these muscles, a single origin-to-insertion path does not encompass their possible actions. One way to solve this is to split the muscle into muscle parts, referred to here as muscle elements. Each muscle element represents a region of the muscle (e.g. the superior portion of adductor magnus) and can have its own set of muscle parameters.
Most modelling systems allow for additional coordinate positions that are used to define the curved path that the physical muscles take between bony attachments (e.g. the path sartorius takes across the anterior thigh as it traverses from anterior superior iliac spine (ASIS) to proximal tibia). Via points, the simplest type of path point, are one or more coordinate positions that are intermediate between the origins and insertion points. The muscle assumes a series of straight line paths between the origin, via point(s) and insertion. The location of the via point is fixed relative to the segment's centre of mass and the connection to other segments. With via points, the muscle is required to pass through each specific point at all times. These points allow muscles to pass over other body surfaces (e.g. iliacus wrapping over the pelvic brim). In addition, muscle geometry paths can be altered by the use of ‘wrapping surfaces’ that are designed to account for relative motion between musculoskeletal elements and other soft tissue (e.g. rectus femoris sliding across the soft tissue structures anterior to the hip). These wrapping surfaces are geometric objects (e.g. cylinders, ellipsoids) that are fixed relative to a particular body segment. Muscle paths are not constrained to pass through any particular coordinate point (as with via points), but rather can move freely over the wrapping surface to find the shortest possible path across the surface between origin, insertion or via points that are beyond the wrapping surface. As the proximal and distal segments change orientation relative to each other, the point of contact between the muscle and the adjacent musculoskeletal elements can move (e.g. gastrocnemius relative to the posterodistal femur as the knee flexes). Forces developed from the action of the muscle act on the body segment used to define the geometry (origin, insertion, via points and/or points on wrapping surfaces). Muscles may also be modelled to be more complex in terms of the number of parameters that define force generation. Minimally, a muscle model must include a maximum force value which defines the force when the muscle is maximally activated, but other approximations of muscle force generation are available, e.g. the Hill-type muscle model [64].
The final required model element is a set of virtual motion markers, which is needed to match model motion to motion measured experimentally (figure 1). For motion tracking (described below), the human subject is fitted with a set of markers that are tracked by a camera system. The musculoskeletal model must include one virtual marker for each experimentally measured marker. Each virtual marker is connected to a body segment and defined in terms of that segment's coordinate system. For instance, if researchers affix a motion-tracking marker over the ASIS of a human participant, the model will require a virtual marker connected to the pelvis body segment with coordinates that position it correctly relative to the model ASIS.
2.2. The input data: motion data and ground reaction forces
An inverse dynamics approach to musculoskeletal modelling (described below) requires that the motions (translation and rotation) of the rigid segments (i.e. the skeletal elements) and externally applied forces be known. This process is well established and has been described in detail elsewhere (e.g. [65–67]). To obtain motions, retroreflective (infrared) markers attached to the skin or clothing are frequently used. Each body segment requires a minimum of three markers that are trackable throughout the entire sequence of movements, so this traceability should be checked in all portions of the data collection volume. The three-dimensional ground reaction force (GRF) for a complete stance is also necessary. GRFs are typically acquired through force plates embedded either in the floor or in a treadmill.
2.3. The simulation
2.3.1. Scaling and marker optimization
Model scaling and virtual marker position optimization are the processes of determining the model segment dimensions and virtual marker positions so that the virtual model matches the human subject and experimental protocol (motion marker placement on the human subject) as well as possible [68]. The degree of matching between the model and experimental data is evaluated using a global error metric between experimental and virtual motion marker positions. During the inverse kinematic analysis (described below), musculoskeletal modelling software finds the position of the model, including overall location in space and individual joint positions, that minimizes the sum of squared deviations between the virtual marker set and experimental motion marker data. The goal of marker optimization and body scaling is to minimize this error term.
The human models that are included with many modelling software packages are generic human models. For analysis, the model must match the subject for whom the input data are known in terms of mass, limb parameters (e.g. distances between joint centres and segment moments of inertia) and virtual marker positions [69,70]. Marker position optimization and scaling are intertwined processes because marker positions can, in part, determine segment dimensions. Initial scaling of the generic virtual model achieves a loose match between the overall bodily dimensions of the subject and those of the virtual model. If a virtual model from a repository is used (rather than creating a virtual model from scratch), then the model will have been developed using the parameters of the individual whose data were available (i.e. a generic model). This generic model will almost certainly not represent the dimensions of the subject on whom motion-tracking data and GRFs were collected, so the generic model needs to be scaled [69]. The simplest scaling approach is uniform scaling of lengths and volumes. Total body mass is the most basic dimension needed, and the ratio of the body mass of the subject to that of the generic model is typically used to adjust segment masses and other parameters that are volume-dependent. The ratio of subject to generic model stature, a length scale, is used as an initial estimation of segment dimensions and other dimensions that are length-dependent. Other anthropometric data (e.g. distances between landmarks) can be used instead of, or in addition to, stature. Other scaling methods are also available, including fully manual scaling, using body fat to adjust mass distribution among the segments, and hybrid scaling (e.g. manually scaling some parameters and algorithmically scaling others). Furthermore, depending on the sophistication of the software used to complete the mechanical analysis, these initial scalings may be adjusted during marker optimization.
Different software packages conduct marker optimization and secondary scaling of segments using different algorithms and require varying levels of interaction/input from the user. At one end of the user-input spectrum, some systems require that the user optimizes the locations of the markers on the model ‘by hand’, moving them relative to their parent segment and then interrogating error values associated with each marker during motion tracking. This process is tedious and may require significant effort before optimal marker position is achieved. At the other end of the spectrum, some software packages change both segment dimensions and virtual motion-tracking markers automatically in a single cohesive optimization process that requires little user interaction. In these frameworks, the user can permit or restrict changes in virtual marker position in none, some or all of the coordinate directions. For instance, motion-tracking markers associated with easily palpable bony landmarks (e.g. medial malleolus) are also easy to locate on the virtual model (referencing the virtual skeleton) and may be considered ‘fixed’ in their location relative to their body segment (the shank). Thus, during the optimization process, such a fixed marker would not change its position relative to the segment, but would influence the determination of shank length. Other markers (e.g. marker plates attached mid-segment) may be allowed to optimize their position in all three coordinate directions, and then contribute little to segment length determination. Finally, some markers may be well defined in two coordinate directions, but less so in the third. For instance, the superoinferior and mediolateral positions of an ASIS marker may be well defined and constrained from moving relative to the pelvis during the optimization process. The anteroposterior position of the virtual ASIS marker relative to the virtual pelvis may be more difficult to estimate as it reflects tissue thickness over this bony landmark (which is hard to measure on living subjects in the absence of internal imaging modalities). Thus, the user might fix the ASIS marker in two coordinate directions, but allow its position to optimize in the third. The fixed mediolateral and superoinferior position on the two ASISs would then influence the scaling of pelvic width.
2.3.2. Inverse kinematics/motion tracking
To carry out inverse dynamics and, in turn, estimate internal forces and moments, segment linear and rotational accelerations must be obtained. Most musculoskeletal modelling software (e.g. Anybody Modeling System™, OpenSim) uses a kinematic analysis procedure of over-determined biomechanical systems, which is formulated as a weighted least-squares optimization problem that matches experimental and model marker positions (e.g. [47,71]).
2.3.3. Inverse dynamics
While statics is the branch of mechanics that deals with stationary objects, dynamics is the branch in which motions are the focus. In statics, an equilibrium that greatly facilitates examination of the system exists where the environment (e.g. the substrate) and the structure (e.g. the human body) are balanced. In standing, the body pushes through the feet and against the ground exactly as much as the ground pushes against the body through the feet. In dynamics, the balance is time-dependent and the solution to the system is through the equations of motion, which relate forces to motions. When the forces that produce the motions are fully known, the problem is characterized as forward dynamics, while when motions are known it is described as inverse dynamics.
In human bodies, the measurement of motions (via marker tracking or other means) is much less invasive and more feasible than measuring the muscle forces that produce the motions, so currently most analyses that depend on the musculoskeletal forces generated in human gait are typically treated as inverse dynamic problems (e.g. [72–74]). Further complicating the solution in multi-body problems (i.e. those with multiple, linked rigid bodies) is that the system can be interrogated from an external or internal perspective. The external view includes the action of the GRF to produce the movement of the body along a path, while the internal one includes both the joint reactions and the action of muscle forces in creating changes in joint angles. In mechanics, the relationship between the number of unknown variables and the number of known variables determines the difficulty of the solution. Fortunately, the internal and external perspectives can be exploited to develop solutions. For instance, in the single-stance foot, the GRF is the external force and produces the moment acting on the foot segment and reacted to by the ankle joint reaction force and moment (which are an internal force and moment). Once the ankle joint reaction force and moment are known, they can be combined with the motion of the shank to solve for the knee joint reaction and moment. The inverse dynamic solution allows for this time-dependent chain solution to determine the joint reactions.
2.3.4. Muscle redundancy solution
A system where the number of variables with unknown value equals the number of independent equations that fully define the system is determinant and a unique solution exists. For systems where there are more variables with unknown values than equations—called indeterminate—assumptions are made to provide more equations [75,76]. Many biological problems, including developing muscle forces, are indeterminate because there are more muscles capable of acting across a joint than necessary from a purely mechanical perspective, i.e. some muscles are (mechanically) redundant. If the desired result from an analysis that uses inverse dynamics is the development of muscle forces, then the joint moments determined from the inverse dynamic solution need to be apportioned to individual muscle forces. This is not a trivial undertaking, because the degree of redundancy is often quite large, so no general consensus on the most accurate criterion on which to base an apportionment algorithm currently exists [77].
Musculoskeletal modelling solutions typically include an algorithm that uses muscle parameters to apportion the joint moments to particular muscles (as muscle forces). Early solutions used PCSA to solve the muscle redundancy problem [75]: the apportionment was made based on the ratio of the PCSA of a particular muscle relative to the total PCSA of all muscles capable of acting (e.g. [61]). While effective from an algorithmic perspective, PCSA-based apportionment does not reflect the natural solution, because cross-sectional area approaches assume equal activation of all possible muscles.
Another algorithm for solving the muscle redundancy problem attempts to minimize the sum of muscle activations raised to a power [76]. The activation represents a proportion of total maximum strength the muscle can generate given the current state of the model (i.e. activation is generally bound between 0 and 1). If the muscle activation exponent is 1, i.e. a linear combination, force is allocated to the muscle with the greatest moment arm until it reaches its maximum allowed value, and then other muscles are recruited. Exponents greater than 1 will distribute forces more evenly across muscles available to produce the moment at a particular joint. The greater the power, the more even will be the distribution (PCSA-based activations are essentially an exponent of infinity). Other criteria for apportionment of the joint moments have also been advanced, such as electromyography (e.g. [78,79]), but activation-based methods remain the standard [76]. Until in vivo muscle forces can be reliably measured, all solutions to the muscle redundancy problem should be treated as hypotheses about how joints and muscles produce motion.
2.4. The output
A valuable aspect of musculoskeletal modelling (or of any model) is that variables of interest can be interrogated anywhere within the model. Standard output from a musculoskeletal model includes linear and angular displacements, velocities and accelerations for all body segments. Similarly, joint positions, velocities and accelerations are also tracked and can be exported for analysis. Muscle length, activation, moment arm and force are all possible outputs, as are joint reaction forces. Additionally, users can track the position of muscle origins/via points/insertions throughout the modelled motion. Conveniently, joint reaction forces can be output in the global coordinate system or the coordinate system of the relevant body segments.
Human locomotion requires the repetition of cycles of limb and body motions that can be reduced to a stride. For walking, a step is defined from an event (e.g. initial contact) on one foot to that same event on the contralateral foot and a stride is two contiguous steps. A stride can vary in duration, even at a constant velocity, so creating an average stride or evaluating the variability of stride parameters requires normalization of the strides to some consistent duration. The typical choice is per cent of stride cycle, where the stride duration is described in terms of equal portions of the entire stride (e.g. [80,81]). Evaluating the variability of parameters of interest (e.g. muscle force, joint reactions) can be an important step in validating the calibration of a solution. While some differences among trials of the same individual are expected, an evaluation of the variability in light of the question that the model is designed to answer should occur.
2.5. Human variability
A crucial advantage of musculoskeletal modelling is the opportunity to include aspects of morphological variation, including the production of subject-specific models. Model scaling and motion tracking incorporate some aspects of human variability, including differences in body mass, limb lengths/dimensions and kinematics, but other aspects are less obviously accessible for modification in musculoskeletal modelling software. These less accessible aspects, such as variation in bone or muscle morphology, are also important to consider. Quantifying human morphological (musculoskeletal) variation is at the heart of biological anthropology. For example, human (and hominin) pelvic morphology has a profound influence on bipedal capabilities [14,80], and has been the focus of intense scrutiny [41,81,82]. Humans are known to be sexually dimorphic in pelvic size and shape [83–85]. Furthermore, human pelvic morphology is also known to vary by population [86–88]. Sources of variation can be incorporated into musculoskeletal modelling to understand pathology [31] and differences in joint reaction forces [35]. For example, marker-informed parameter optimization and secondary scaling (described above) can change the width of the pelvis if experimental and virtual markers are included on bony landmarks of the pelvis, such as the ASISs. If the subject is wider than the generic model, the parameterization can widen the pelvis so that the virtual and experimental markers are more closely aligned. Scaling can also be accomplished manually either by applying factors to the three segment coordinate directions or by redefining the locations of specific anatomical landmarks associated with the segment. The former method moves all points proportionately while the latter allows for localized warping.
In addition to skeletal variation, aspects of muscle morphology that vary among humans may also be critical for accurate musculoskeletal modelling. For instance, variability in muscle parameters has been of concern for some time [89]. While repository models represent a particular person's morphology, people differ in important musculoskeletal modelling parameters such as muscle volume, optimal fibre length and pennation angle (table 1).
Table 1. Comparison of muscle parameters between the standard MoCap model in AMMR™ v. 2.3 (AnyBody Technology, Aalborg, Denmark; from Klein Horsman et al. [90] and Carbone et al. [34]) and Charles et al. [91]. Muscle names are per the MoCap AMMR v. 2.3. Where multiple muscle regions are present in the MoCap model, the values for the muscle volumes are summed while the optimal fibre lengths and pennation angles are averaged.Collapsemuscleoptimal fibre length (cm)pennation angle (°)muscle volume (ml)MoCapCharles—averageMoCapCharles—averageMoCapCharles—average
adductor brevis | 10.38 | 7.6 | 0.0 | 11 | 108.90 | 93 |
adductor longus | 10.57 | 11.0 | 0.0 | 12 | 159.56 | 159 |
adductor magnus | 9.97 | 23.1 | 0.0 | 12 | 559.25 | 567 |
biceps femoris caput breve | 9.14 | 10.9 | 0.0 | 9 | 107.95 | 92 |
biceps femoris caput longum | 8.54 | 20.4 | 29.9 | 11 | 232.01 | 194 |
extensor digitorum longus | 5.99 | 13.8 | 8.3 | 7 | 32.29 | 76 |
extensor hallucis longus | 5.98 | 10.6 | 14.4 | 7 | 36.27 | 21 |
flexor digitorum longus | 3.84 | 28.4 | 25.28 | |||
flexor hallucis longus | 5.20 | 30.1 | 161.50 | |||
gastrocnemius lateralis | 5.69 | 12.2 | 25.4 | 9 | 136.36 | 128 |
gastrocnemius medialis | 6.01 | 9.7 | 10.8 | 10 | 263.26 | 230 |
gemellus | 3.43 | 0.0 | 28.40 | |||
gluteus maximus | 13.57 | 0.0 | 936.55 | |||
gluteus medius | 6.76 | 5.3 | 674.71 | |||
gluteus minimus | 3.28 | 0.0 | 82.59 | |||
gracilis | 18.11 | 17.3 | 0.0 | 6 | 87.97 | 91 |
iliacus | 8.16 | 8.8 | 203.13 | |||
obturator externus inferior | 6.88 | 0.0 | 37.88 | |||
obturator externus superior | 2.77 | 0.0 | 68.18 | |||
obturator internus | 2.05 | 0.0 | 52.08 | |||
pectineus | 9.50 | 0.0 | 64.58 | |||
peroneus brevis | 4.54 | 23.1 | 86.09 | |||
peroneus longus | 5.08 | 15.8 | 121.52 | |||
peroneus tertius | 4.27 | 19.1 | 26.52 | |||
piriformis | 3.88 | 0.0 | 31.25 | |||
plantaris | 4.77 | 0.0 | 11.36 | |||
popliteus | 2.40 | 7.4 | 0.0 | 8 | 25.57 | 15 |
psoas major | 9.92 | 13.4 | 193.18 | |||
psoas minor | 5.91 | 0.0 | 6.63 | |||
quadratus femoris | 3.37 | 0.0 | 49.24 | |||
rectus femoris | 7.84 | 14.2 | 21.9 | 8 | 226.33 | 249 |
sartorius | 34.71 | 40.8 | 0.0 | 205.49 | 143 | |
semimembranosus | 8.09 | 15.8 | 25.0 | 12 | 138.26 | 244 |
semitendinosus | 14.16 | 18.3 | 0.0 | 8 | 208.33 | 186 |
soleus | 4.40 | 14.6 | 61.6 | 12 | 792.91 | 461 |
tensor fasciae latae | 9.49 | 0.0 | 83.33 | |||
tibialis anterior | 6.83 | 13.7 | 9.6 | 7 | 181.49 | 129 |
tibialis posterior | 3.78 | 34.2 | 163.36 | |||
vastus intermedius | 7.68 | 18.1 | 11.8 | 11 | 292.61 | 521 |
vastus lateralis | 6.69 | 19.6 | 0.0 | 15 | 583.33 | 606 |
vastus medialis | 7.16 | 15.9 | 0.0 | 14 | 454.27 | 415 |
Several possibilities for future work that explores the impact of morphological variation on muscle and joint forces are worth noting. For instance, while some anthropological work has elucidated the interactions of morphology and walking conditions (e.g. velocity, changes of direction, burdens and gradients) on energy expenditure (e.g. [92–94]), much less is known about how morphology and walking condition factors interact to affect joint and muscle forces (although see [31] for an example). People vary morphologically and walking under varied conditions is an activity of daily living, so this area is an important one for many disciplines to pursue including biological anthropology, biomedical engineering and product design. Perhaps more importantly, collecting the requisite data to inform a musculoskeletal model is time consuming, invasive and dependent on (at this time) fragile laboratory equipment. Consequently, acquiring samples that represent the breadth of human morphological variation is unlikely to occur in the near future. Once validated, though, a musculoskeletal model can be modified to assess the impact of these known morphological differences on muscle and joint forces. Extending that line of research further, a longer-term goal would be use musculoskeletal modelling to understand the influence of the morphological differences between modern humans and extinct hominins on muscle and joint forces, an area that has sought biomechanical solutions without effective software support for decades [12,14,61,95,96].
2.6. Current study
Here, we demonstrate the process and output of one human subject during walking using AnyBody Technology's musculoskeletal modelling software (Anybody Modeling System™; AnyBody Technology, Aalborg, Denmark). We present several of the possible outputs that could be of interest to biological anthropologists. This includes joint kinematics, joint reaction forces and a selection of lower limb muscle force magnitudes. We make some recommendations for future researchers to help facilitate their use of software and output. Additionally, by presenting a set of possible outputs, we demonstrate the motivation for undertaking musculoskeletal modelling and that results may have a direct bearing on anthropological questions, or as intermediate steps to other analyses (e.g. finite-element analysis). Additionally, while engineering has much to offer biological anthropologists, the latter have an appreciation for human variation which may enhance future generations of musculoskeletal models.
3. Material and methods
3.1. Participant
One healthy participant was used for this demonstration, but he is part of a larger study aimed at examining human walking. The participant was 36 years old (male; body mass, 74.3 kg; stature, 172 cm) and reported being free from lower-limb injuries. The University of Washington's Institutional Review Board approved all aspects of this study (IRB no. STUDY00001125).
3.2. Experimental protocol
Kinetic and kinematic data were measured using a four-force plate (Kistler, Switzerland), 10-camera motion capture system (Qualisys, Sweden) in the Amplifying Movement & Performance (AMP) laboratory of the University of Washington. Thirty infrared-retroreflective 5 mm hemispherical markers were affixed to anatomical landmarks and used to track motion through the laboratory (figure 1; electronic supplementary material). The participant walked unshod the length (10 m) of the gait laboratory five times at his self-selected preferred walking speed (1.15 m s−1). The participant walked in a straight path with his eyes directed at a point on the far wall of the laboratory and contacted the surface of each force plate with a single foot. Trials with multiple feet on a single force plate or ones where the stance foot crossed plates were immediately redone. The participant was allowed several attempts prior to data collection to become familiar with the protocol. All trials exhibited similar kinematics and kinetics (data not shown), so we selected one trial to use in this demonstration.
Marker data and GRFs were collected at 120 Hz and 1200 Hz, respectively, and then were filtered at 5 Hz and 10 Hz, respectively, with a fourth-order, low-pass zero-lag Butterworth filter. Calibration of the system yielded a limitation in its fidelity for marker data of 1 mm and force data of ±2.5 N for the direction of travel (X), ±5 N side (Y) and ±25 N vertical (Z). All data were exported from the Qualisys Track Manager software in C3D file format, which can be read directly by the modelling software (provided in electronic supplementary material).
3.3. Musculoskeletal model
We used the ‘MoCap’ model from the AnyBody Managed Model Repository™ (AMMR v. 2.3; AnyBody Technology, Aalborg, Denmark; http://www.anybodytech.com/software/model-repository-ammr), which is a multi-trial, full-body, motion-capture-driven human gait model, and the commercially available AnyBody Modelling System™ (v. 7.3; AnyBody Technology, Aalborg, Denmark) to calculate forces in the lower limb [34,97]. The MoCap model that we used for this demonstration has been used to assess human gait in numerous applications (e.g. [98,99]). The MoCap model is assembled from a trunk and left and right lower limb components. The trunk includes a lumbar region, trunk, neck and head. The lower limbs consist of the following segments: thigh, patella, shank, talus and foot. Each lower limb has six total DOFs, including all three rotations at the hip and one each at the knee (flexion/extension), ankle (plantar flexion/dorsiflexion) and subtalar (inversion/eversion) joint. The pelvis (relative to the ground) has six DOFs, three translational and three rotational. Forty-one anatomically defined, bilateral muscles (table 1) are represented by 169 muscle elements in each model lower limb (e.g. gluteus medius is composed of 12 separate muscle element actuators) [34]. The muscle elements in this model generate force as a function of their activation and a strength parameter which is appropriate for human walking [97,100]. We defined a virtual motion-tracking marker for each of the experimentally measured motion markers (electronic supplementary material). The polynomial algorithm that was employed to solve the muscle redundancy problem apportioned force to muscles using an objective function with power = 3. (For more information regarding the model, see §5.3.)
3.4. Musculoskeletal simulation
The first step of the simulation was to scale the model to match the dimensions of the subject while simultaneously optimizing the marker coordinate and joint rotation axes' locations [68]. After initial scaling, the values of these parameters (i.e. segment dimensions, marker locations, rotation axes' locations) were determined using the marker motion data from a trial. Values for all parameters were optimized to minimize the sum of square distances between model marker positions and experimental marker positions. Following convergence of parameter optimization, we conducted an inverse kinematics analysis. Finally, we conducted the inverse dynamic analysis, which includes apportioning muscle forces using the algorithm that minimizes the cost function as a sum of muscle activation values raised to the third power. This value has been shown to allocate muscle activations in a way that reflects experimental data [77,100,101].
3.5. Output variables and data processing
We output several kinematic and kinetic variables, including joint kinematics, filtered GRF components and joint reaction forces (electronic supplementary material). Additionally, we output the force magnitude for a selection of the muscle elements for the left lower limb during an entire stride cycle of walking. Muscle element forces within anatomically defined muscles were summed for presentation (e.g. 12 elements that represent gluteus medius). All time-related curves were resampled to 1% increments of the stride cycle.
4. Results
GRF component curves across the gait cycle are presented in figure 2. As is typical for human walking, the vertical component has two peaks (early stance peak, 719 N; late stance peak, 778 N) on either side of a valley (633 N) which occurs around the middle of stance phase. The average of the two peaks is 103% of body weight (729 N), while the valley represents 87% of body weight.
Figure 2. Filtered ground reaction force components. Grey vertical line represents moment of toe-off.
Pelvic rotation and obliquity and hip, knee and ankle kinematics are provided in figure 3. All kinematics are of the left lower limb throughout a single stride. Pelvic motion is described in terms of the right hip (previous stance limb) relative to the left hip (current stance limb) [102]. Hip, knee and ankle joint reaction forces, which include contributions from muscle forces, are provided in figure 4. In each joint, the vertical joint reaction force is greater in magnitude than the anteroposterior and mediolateral components. All three joints show a large peak during the second half of the stance phase, and this late peak is significantly larger in the ankle. Muscle force profiles for 18 anatomically defined muscles across the stride cycle are provided in figure 5. During the stance phase of walking, gastrocnemius, soleus and gluteus medius generate the largest forces. All output data are provided in the electronic supplementary material.
Figure 3. Angular kinematics for the pelvic segment and hip, knee and ankle joints. Vertical grey line represents toe-off.
Figure 4. Joint reaction forces expressed in the global coordinate system. Vertical grey line represents toe-off.
Figure 5. Muscle forces for 18 anatomically defined muscles.
5. Discussion
5.1. Model output
The subject that we used to demonstrate the musculoskeletal modelling approach to modelling walking is an adult male without gait pathology and within species-typical anthropometric characteristics. Consequently, the model-derived motions should accord with those determined experimentally and extensively described in the last 40 years of gait research (e.g. [66,103–105]). The model results presented above do meet this requirement, as discussed below in detail. We focus on the model-derived kinematics, because these kinematics are the foundation for all subsequent analyses in musculoskeletal models. Any analysis should begin by verifying that model-derived variables align with those obtained empirically, and the cause of any discrepancies should be understood or noted as a potential limitation. We belabour the details in the following paragraphs in order to demonstrate the basic approach to model validation.
In our subject at heel strike, the pelvis is rotated approximately 10° in the horizontal plane such that the right hip is posterior to the left hip (e.g. the current stance limb). The pelvis then swings through a neutral position and, while leading up to the next right heel strike, into an orientation with the right hip anterior to left hip (approx. 10°). The pelvis drops in the frontal plane with the right hip inferior to the left (stance) hip soon after heel strike as the right lower limb moves into the swing phase with pelvic obliquity ranging between ± 4°. The pelvis remains in this position until the subsequent right heel strike, after which the pelvis approaches the neutral position. Model-derived pelvic rotation and obliquity are similar to empirically determined values (e.g. [103]).
In the sagittal plane, the left hip begins the stance phase in flexion of approximately 17° and moves through a neutral position and into extension of approximately 17°. Winter [66] describes a less symmetrical flexion–extension cycle at natural cadence (mean values of approx. 20° of flexion and 11° of extension), but the model-derived hip rotations are well within one standard deviation of his mean values. During most of the stance phase, the left thigh is adducted less than 6° in accord with others (e.g. [103]).
While the pelvic and hip joints have three rotational DOFs, the knee and the ankle only rotate in the sagittal plane. The left knee flexes slightly (10°) early in the stance phase, but maintains a near fully extended position through midstance, after which it begins to flex in preparation for the swing phase where flexion exceeds 60°. After the initial heel strike and as the forefoot contacts the substrate, the ankle plantar flexes slightly (less than 5°) and then assumes a dorsiflexed position for most of stance. Rapid plantar flexion occurs in late stance as the body is propelled forwards, reaching a plantar flexion value of greater than 10°. The ankle returns to a dorsiflexed position for swing. These values for both knee and ankle angles are typical (e.g. [66,103,105]).
Consequently, the kinematics of the model match those that we expect empirically. While this check is not sufficient to ensure that the internal (muscle and joint reaction forces) are reasonable, it is a necessary condition and validation for any model. Incorrect kinematics guarantee incorrect or suspect muscle forces.
Joint reaction forces have been empirically measured only with instrumented joint implants. By definition, the surgery creates conditions that limit applicability of the results, but no other method currently exists to measure force. Pedersen et al. [45] measured hip joint contact forces in a 72-year-old man (63.2 kg) with a hip replacement and found that the vertical force was 1750 N, considerably less than the almost 3500 N we found with the musculoskeletal model. Important differences between this patient and our subject exist, including age (72 versus 36 years), body mass (62.3 versus 74.3 kg), health status (58 days post hip replacement surgery versus healthy) and walking speed (0.89 versus 1.15 m s−1). While hip and knee joint reaction forces presented here are larger than in vivo values using instrumented prostheses [77,106,107], our model-derived joint forces are, however, similar to other models of healthy adults [35,107].
Although it is difficult to compare muscle forces with experimental data, it is possible to consider the pattern of activation with reference to the muscle's action. The muscles that generate the highest force in the model are gluteus medius, soleus and gastrocnemius. Gluteus medius plays a critical role in bipedal walking by preventing pelvic drop during single stance [103,105]. The gluteus medius muscle profile from the model is double peaked, and the timing of the two peaks is coincident with the two peaks in the vertical GRF, indicating the muscle's support of the pelvis during single stance. Gastrocnemius and soleus plantar flex the ankle; when the foot is in contact with the ground and the hip is extended, plantar flexing the ankle moves the body anteriorly [103,105]. The gastrocnemius and soleus muscles show the greatest force generation during the propulsive portion (second half) of the stance phase, in keeping with our expectation from their function.
Other muscles also display moderately high forces. Tibialis anterior, which dorsiflexes the ankle, shows muscle force generation immediately after heel strike when the body begins to move over the foot, and then is mostly inactive until the second half of the swing phase. The quadriceps demonstrate a force peak early in stance phase, and again at the stance–swing transition, potentially relating to their role in extending and stabilizing the knee [105]. The hamstring muscles show the greatest peak at the end of swing phase and into early stance phase. Activation of these muscles at the end of swing phase is associated with decelerating the limb prior to heel strike. Early in stance the GRF is directed anterior to the hip and knee, and the hamstring muscles are active to prevent further flexion (hip) and extension (knee) of these joints. The hip flexors are active late in the stance phase, and may be active to counteract hamstring hip extension moment as the hamstrings flex the knee in preparation for swing phase. The adductors show initial peak in stance (adductor magnus) and late stance (adductor longus), but generally smaller than forces in other muscles.
5.2. Input data choices
The model-derived motions and joint reaction and muscle forces align with our expectations from the experimental data in the literature. We achieved these results through careful attention to detail when we collected the data we used as inputs to the musculoskeletal model. We provide below some brief suggestions from our experience. Because motions are derived from the change in position of these markers, marker placement on the segment should be chosen to define fully the motion of the segment and to reduce marker positional noise as much as possible. Positional noise arises from the inherent error of measurement of the marker's position in space (system error) and from movement of the surface (skin or clothing), where the markers are placed relative to the underlying rigid segment (e.g. [108–110]). System error is reduced through calibration; for our system, it is less than 1 mm for any marker. Movement between the external (skin) surface and the bone is more complicated to reduce, but may not be clinically relevant [110].
Marker position is a compromise with soft tissue coverage of underlying bone, landmark palpability (to aid placement) and distance between markers (such that marker positional error is much less than the distance between markers) all being important. The shank and foot have palpable landmarks that are covered by minimal soft tissue (e.g. the malleoli), but the thigh (e.g. greater trochanter) and, especially, the pelvis are more challenging segments on which to affix markers. This problem is exacerbated when the subjects of interest have more soft tissue (muscle or fat) overlying the bony landmark, and considerable soft tissue coverage is common in the hip–pelvis region. One approach to ameliorate some of this difficulty is to use a marker plate which has four (or more) markers mounted on a hard plastic plate that is then strapped to the segment. If a plate can be affixed firmly to the segment, as it can be to the thigh, the marker plate eliminates between-marker error. Unfortunately, marker plates are difficult to affix firmly to the pelvis. Marker placement issues, consequently, are most difficult to overcome to track the motion of the pelvis. Yet, for musculoskeletal models of human locomotion, pelvic motions are critical components of system performance. We achieved the motions shown in figure 3 through extensive protocol development and testing.
Another critical input to the simulation are the GRFs. It is possible to conduct an inverse dynamics solution of the single support phase of walking with a single force plate (e.g. [111]), but extending the analysis to the double support portions of the stride is more difficult with one force plate because it would require estimation of the unknown contribution of the other stance foot. Consequently, any research that requires information about double support should optimally be conducted with multiple force plates. The placement of the force plates relative to each other should be appropriate to the population and gait. Stride lengths for adults or running gaits are longer than for children or walking, and plates should be positioned to prevent subjects ‘targeting’ foot placement or large numbers of repeated trials. Data analysis is less complicated if stance occurs on one plate, but the data from two immediately adjacent plates can be used as long as there is no gap between them [112]. Multiple feet contacting a single plate is also possible to resolve with additional information (such as localized pressure sensors) or predictive methods (e.g. [113,114]). All of these remediation techniques are time-intensive. The laboratory in which the presented data were collected has four floor-embedded force plates that are positioned to allow normal stride lengths for adults when walking. The subject and trial that we used for this demonstration is part of a much larger study [115]. Nonetheless, we assessed each foot placement on all four plates for each trial before moving on to the next trial. Such data checking is a tedious effort at the time of data collection, but well worth the effort to ensure all data are suitable for future analyses. Using four force plates, we had data frames before and after the data that we used in the analysis. This allowed us to discard the initiation and termination frames of the simulation, which are frequently unusable. Consequently, researchers are advised to locate force plates well and test the placement prior to data collection.
Finally, researchers may consider collecting complete anthropometrics of subjects and using these in the initial scaling as the starting points for the parameterization. People vary considerably in terms of their proportions [35,37], so standard models, such as the model we used here, represent a morphology that may not match the dimensions of any particular subject. If your model starts with segment lengths that are substantially different from the segment lengths of the model, the parameterization algorithm may not be able to converge on a suitable solution. Two failures can occur: the parameterization can fail to converge, which produces an obvious error message, or the parameterization converges, but, in order to converge, the simulation requires changes to the morphology or motions that are inconsistent with the experimental data or subject dimensions. This latter failure can be insidious. An absolutely critical step in model validation is to check the motions and shapes of the segments after parameterization. A poorly parameterized model is incorrect and will not produce acceptable motion or joint and muscle forces.
5.3. Modelling choices
In order to produce this demonstration, we made a number of choices with regard to the model which bear discussion. The most important of these was the choice to use a standard, repository model, in this case the MoCap model of AMMR v. 2.3. By using the MoCap model we accepted many assumptions and compromises, but gained the experience and efforts of many dedicated musculoskeletal modellers. Nonetheless, the responsibility of assessing the ability of any model to answer the question of interest rests with the researchers using the model. Our intention was to demonstrate some of the decisions and validations needed in musculoskeletal modelling of human walking and our choices reflect that.
The first of these choices is the theory on which the muscle actuation rests. The muscle models used here generate force as a function of maximum force and activation. More sophisticated muscle models, such as the Hill muscle model, are available. The Hill muscle models are thought to provide the most complete model of muscle behaviour [116], but the Hill model requires several parameters for its full definition, including: maximum isometric force, tendon slack length, optimal fibre length, pennation angle at optimal fibre length, activation and deactivation times, as well as maximum contraction velocity. The number of parameters that must be defined for a Hill muscle model, as well as the sensitivity of muscle behaviour to parameter values, means that considerable effort must be expended to parameterize and optimize Hill-type models. As with all models, poorly informed inputs (i.e. muscle parameters) lead to unstable results. This means that for certain activities it may be advisable to use a simpler, robust model rather than a more complex, unreliable one. For human activities during which muscles operate at slow contraction velocities and near their optimal fibre length (e.g. human walking), simple muscle models that depend only on geometry and a maximum force parameter are preferable [100]. Consequently, we accepted the use of the simple muscle model employed in the MoCap gait model.
Another choice was the selection of the scaling method. Understanding which scaling method to use (e.g. mass–length–fat), requires understanding the underlying assumptions of the method (e.g. distribution of mass to segments), to which model parameters the scaling applies (e.g. which segment length scale to which muscle length) and how the model uses the scaling to define parameters (e.g. are all parameters scaled or just a subset?). We made many other choices by not changing components that are inherent to the generic MoCap model, including acceptance of the form of the algorithm used to solve the muscle redundancy problem. Such choices are often embedded in standardized models, and the user does not have to modify any of them in order for the model to run to completion. This makes ensuring that the question of interest can be addressed by the model become an intentional assessment. For instance, the standard MoCap model has only one rotational DOF for the knee, which works well to assess overall walking parameters but is less well suited to a detailed study of the knee. Furthermore, standardized models are improved as more experience and data become available (e.g. [97]). For instance, the MoCap model's leg (TLEM v. 2.1) has been modified to reflect muscle paths more accurately by the inclusion of wrapping surfaces.
We also relied on the standardized definition of available outcome variables such as muscle forces and joint motions. Given that most software allows for forces and motions to be expressed in different ways (at a minimum, in different coordinate systems) it is imperative to understand the context of the reported variable. Is an output reported in the global or segment coordinate system? Is a joint force an internal force or a reaction? In what direction does the muscle force act? As with all big data, extensive data exploration and a complete understanding of the theories (mechanics), procedures (inverse dynamics) and systems (human anatomy) is necessary to draw conclusions from a musculoskeletal modelling simulation.
5.4. The role of anthropology in musculoskeletal modelling
Finally, it is important to make clear that in this demonstration we have used a walking trial of one subject and focused on the subject's left side. The kinematics and GRFs produced by this individual's right side and by other walking trials demonstrate a consistent pattern, but were not numerically identical to those reproduced here (electronic supplementary material). Consequently, the muscle forces and joint reactions are also somewhat different. The degree of intraindividual variability in gait parameters is influenced by such characteristics as age and health, but, even for this healthy adult male, peak forces within a muscle can vary. In muscles that exert relatively small forces (less than 500 N), step-to-step values can vary by more than a factor of 2 (electronic supplementary material). Muscles which show greater force generations (e.g. gastrocnemius; greater than 1000 N), step-to-step variability is lower (a factor of approx. 1–1.4; electronic supplementary material). This serves to underscore the importance of collecting and analysing multiple steps when the intent is to report patterns and forces (not the intent here).
While scaling addresses some aspects of variability and improves model predictions [35,69], scaling itself does not solve all the problems. Interindividual variability in such important variables as muscle parameters has been appreciated for over a decade [89], but the consequences of such variability for musculoskeletal modelling remains unclear. For instance, muscle strength as implemented in the solution that we use for this report depends on three factors: muscle volume, optimal fibre length and pennation angle. The MoCap model, which we used in our demonstration, uses a set of values derived from several sources but is critically informed by Klein Horsman et al. [90] and Carbone et al. [34]. These data are from measurements of cadavers who were 77 and 85 years old at the time of death, respectively. Cadaveric data are widely used in musculoskeletal models (e.g. [53,99,117]), but muscle parameters from cadavers may not represent those from individuals who are not at imminent risk of death. Table 1 includes values from Klein Horsman et al. [34] and Carbone et al. [34] and from a recent analysis of 10 healthy individuals who were evaluated with diffusion tensor imaging [91]. Gastrocnemius and soleus, the principal drivers of forward propulsion, vary in all three factors. Importantly, the complexity of gait mechanics that makes musculoskeletal modelling an attractive solution also makes understanding the effect of this variation impossible without simulating the difference. Larger volume, smaller pennation angle and shorter optimal fibre lengths all act to increase muscle strength, thereby increasing the allocation of force to it by the algorithm. But, the algorithm uses a muscle's strength relative to all other muscles' strengths, making it impossible to predict the impact of muscle parameter variation on muscle forces or joint reactions. Variation in muscle parameters is but one of the many ways that people vary in their musculoskeletal anatomy. Leveraging the expertise of biological anthropologists to assess the impact of human biological diversity on musculoskeletal modelling can only be beneficial for the entire community of researchers.
6. Conclusion
Although motivated by different specific goals, biological anthropologists and engineers (principally biomedical) both pursue questions that revolve around understanding the internal forces of the primate musculoskeletal system. Engineering questions tend to be species-specific, focusing on living humans. Biological anthropologists are also interested in living humans, but frequently ask questions that include greater taxonomic diversity. Research in both fields can be specific to the details of the individual (e.g. patient-specific orthopaedics, particular hominin fossils) or generic representatives of groups (e.g. sex-specific running shoes, comparisons of primate taxa).
Musculoskeletal modelling is powerful because of the flexibility it provides at every stage of research. Building a new model requires intensive effort, but may be necessary based on the question. Previously validated models from repositories lower barriers to modelling experiments, but researchers should be aware of choices made by model creators. Like all modelling exercises, serious and intentional thought must be dedicated to the multitude of modelling choices. This is not to say that all models need to be complex. If a simpler model adequately describes the phenomenon of interest and answers the question at hand, then a simpler model may be preferable.
Critical to support rigorous simulations, musculoskeletal models include changeable parameters that describe all aspects of the musculoskeletal system. This is a purposeful feature included by engineers that drove musculoskeletal modelling and an attractive feature for biological anthropologists. Model users can change the size and shape of the skeletal elements, dimensions and inertial properties of body segments, DOFs and range of motion of joints, action of passive soft tissue structures, the geometry and force capacity of muscles, as well as numerous other parameters. Generic models from repositories (such as the one used here) can be scaled and reshaped to capture many aspects of human variation (e.g. limb dimensions, mass), and hence are well suited to many questions. More effort can be dedicated to refine aspects of a model that are of particular interest (e.g. change the morphology of the bony pelvis). Modelling flexibility naturally extends to the information that can be derived from simulations. Forces generated by individual muscles (or portions of muscles) as well as internal to joints are standard output, as is the ability to express these quantities in different coordinate systems. Until internal forces of the body can be measured non-invasively, musculoskeletal modelling will remain the pre-eminent approach to these questions.
Future collaborations between engineers and biological anthropologists to create better models of human movement have the potential to expand the understanding of mobility in both fields. Subject-specific models are internally consistent but not generally representative of people, so extrapolating should only be done with this limitation in mind. Developing a model from scratch for a fundamentally different organism (e.g. a chimpanzee) is a significant investment.
Data accessibility
Additional information is available as the electronic supplementary material. Raw data and simulation output are available at the assigned figshare link.
Authors' contributions
A.D.S. and P.A.K. conceived the project. S.G.L. and P.A.K. designed and carried out the data collection. A.D.S. and P.A.K. carried out the simulations and data analysis. A.D.S. and P.A.K. wrote the draft and all authors contributed to the final version of the manuscript.
Competing interests
We declare we have no competing interests.
Funding
P.A.K. and S.G.L. thank S.K. Benirschke, MD, Jerome Debs Endowment Chair in Orthopaedic Traumatology at the University of Washington, for the ongoing support of this research and for our many conversations about foot form and function.
Acknowledgements
The authors would like to thank the participants of this study for their time and patience during the protocol. The authors would also like to thank the reviewers and editors for comments that improved this manuscript.
Footnotes
One contribution of 12 to a theme issue ‘Biological anthroengineering’.
Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5486386.
Review articles
A review of musculoskeletal modelling of human locomotion
,
and
Published:13 August 2021https://doi.org/10.1098/rsfs.2020.0060
Abstract
Locomotion through the environment is important because movement provides access to key resources, including food, shelter and mates. Central to many locomotion-focused questions is the need to understand internal forces, particularly muscle forces and joint reactions. Musculoskeletal modelling, which typically harnesses the power of inverse dynamics, unites experimental data that are collected on living subjects with virtual models of their morphology. The inputs required for producing good musculoskeletal models include body geometry, muscle parameters, motion variables and ground reaction forces. This methodological approach is critically informed by both biological anthropology, with its focus on variation in human form and function, and mechanical engineering, with a focus on the application of Newtonian mechanics to current problems. Here, we demonstrate the application of a musculoskeletal modelling approach to human walking using the data of a single male subject. Furthermore, we discuss the decisions required to build the model, including how to customize the musculoskeletal model, and suggest cautions that both biological anthropologists and engineers who are interested in this topic should consider.
1. Introduction
A principal research focus within biological anthropology is interpreting skeletal variation within the context of behavioural diversity, including variation in diet, disease processes and activity patterns among living and extinct primates, including humans and hominins [1–4]. Locomotor behaviour and musculoskeletal morphology have been a central focus of this area of inquiry because locomotion provides primates with access to key evolutionary resources including food, shelter and mates. Furthermore, bipedalism is considered a defining character of the hominin lineage, cementing its importance within the field [5–7]. The challenge of elucidating form–behaviour relationships are the complexities of the intervening and underlying functions (e.g. joint range of motion) and biomechanics (e.g. joint reaction forces).
Early anthropological research connecting form and function was both typological and qualitative in nature; however, anthropologists have increasingly adopted mechanical engineering approaches to answer questions about primate locomotion generally, and hominin locomotion in particular (see [4] for a historical review). Within anthropology, two major subfields of mechanics have been applied to locomotor systems: Newtonian mechanics and solid continuum mechanics. Newtonian mechanics is concerned with the description of the motion of solid bodies under the influence of a system of forces, while solid continuum mechanics quantifies the behaviour (e.g. deformation) of solid materials when subjected to forces [8,9]. Lovejoy and colleagues [10–13] were early adopters within anthropology of both Newtonian mechanics and solid continuum mechanics [4]. Newtonian mechanics is routinely used by anthropologists to investigate a variety of questions (e.g. [14–19]). Beam theory, an application of solid continuum mechanics, has been used to estimate the structural capacity of long bones and is regularly exploited by anthropologists (see [4] for a review) to reconstruct hominin locomotor adaptations (e.g. [13–17]) as well as modern human mobility patterns (e.g. [18–24]).
The application of Newtonian and continuum mechanics is, of course, a core competency of mechanical and civil structural engineering, and the methods that anthropologists use derive directly from these engineering fields. The techniques and theories that underlie their application were developed, however, to design and analyse human-made products using inert materials. While such structures and machines tend to be faithfully produced based on design specifications, humans are not clones. Rather, each individual is completed from a standard template that has, as an inherent feature, the ability to respond to the environment. This means that whereas human-made products generally have low tolerance for variance, among humans, variation is the rule. Furthermore, the ‘standard’ template itself varies among ancestral lineages, populations and species because of the evolutionary history of each (e.g. [25–27]). The biomedical engineering world has responded to this with the use of patient-specific models because, in patients, a direct relationship exists between a particular morphology (i.e. that of the patient) and the question to be answered (e.g. which implant configuration limits strain on the bone?) with the musculoskeletal model (e.g. [28–35]). Not all questions, however, concern specific individuals. In anthropology, the focus may be comparing species (or populations) to understand responses to selective pressures. In engineering, running shoes, for example, are not designed for particular individuals, who vary in foot function [36], but rather for people generally. And people are not just scaled versions of each other [37]. Thus, an appreciation of human variation, a core competency of biological anthropology, is a critical component of any modelling exercise.
At the heart of many locomotor-focused questions, both engineering and anthropological, is a need to evaluate internal (e.g. muscle, joint contact) forces. Finite-element analysis of skeletal elements [38,39], estimating locomotor costs [40,41] and evaluating task performance [42,43] all require muscle forces to be known or estimated. Measuring muscle and joint contact forces in living animals is, however, exceptionally challenging because it requires surgical intervention to implant the relevant sensors into the body [44,45]. As an additional obstacle, these techniques are normally limited to a single joint or muscle/tendon within a single study (see [44] for examples). Consequently, internal forces must be estimated. Poorly estimated loads are likely to produce misleading results, interpretations and conclusions. Consequently, estimating internal forces (e.g. muscle forces) in humans is of critical interest to both anthropologists and engineers.
1.1. Musculoskeletal modelling: a path forwards
Originally driven by a clinical interest in gait pathologies, musculoskeletal modelling has emerged as the pre-eminent technique to elucidate the internal details of human (and other animal) movement [46]. These models build upon dynamic analyses of linked rigid-segment models ([47] and described below). Rigid-segment models represent the body as a series of solid body segments which are linked together at joints. Musculoskeletal modelling extends rigid-segment models by including detailed models of individual muscles, as well as models of neurological control. In extending the dynamic approach, musculoskeletal modelling solves the muscle redundancy problem based on muscle parameters under a specified minimization criterion [48–51]. The muscle redundancy problem refers to the fact that there are more muscles than necessary to actuate movements at joints [48,49]. In solving for muscle forces, this leads to an indeterminate system of equations in which there are more unknowns (muscle forces) than constraints (equations) [51], thus requiring optimization solutions. Results of musculoskeletal modelling can reveal individual sequences of muscle excitation and activation, as well as estimates of muscle forces and bone-to-bone contact forces.
In a series of papers, Pedersen et al. (see [45]) developed an early version of this type of model and used it to estimate lower limb muscle forces and hip joint contact forces. The hip joint contact forces were validated using an instrumented hip joint replacement in a single elderly subject [45]. Consequently, these data provided an early demonstration of the validity of musculoskeletal modelling as an approach to determining internal forces. Since then musculoskeletal modelling has been used for numerous projects exploring human motion and locomotion [31,43,52–56] as well as that of non-human primates [57–59] and early hominins [60]. Wang et al. [61] applied this technique to the hominin fossil record, but their approach was necessarily limited by the lack of software available to carry out the complex calculations. Fortunately, several software packages (e.g. AnyBody Modeling System™, OpenSim, FreeBody) are now available to carry out musculoskeletal modelling [45,62,63].
Accurate musculoskeletal models demand both understanding of the principles of Newtonian mechanics and an astute appreciation of anatomy and anatomical variation. Thus, successful musculoskeletal modelling is clearly situated at the interface between engineering (i.e. mechanics) and biological anthropology (i.e. human variation). Here, we provide an overview of this effort from an anthroengineering perspective, blending knowledge from both domains with a special emphasis on a task—understanding human walking—that is of equal interest to both groups.
2. Background
2.1. The model
In musculoskeletal modelling used to solve for muscle forces, the body of interest (e.g. the whole human body) is modelled as a series of rigid segments connected by joints and actuated by muscles [34,46] (figure 1). This body moves in a virtual space whose dimensions are defined by a global Cartesian coordinate system. Models can include as few or as many body segments as necessary to answer the question of interest. For walking, models minimally include segments for the pelvis (and trunk/arms/head, often aggregated) as well as left and right feet, shanks and thighs. More complete models often also include segments that represent the patella and talus, as well as a head–arm–trunk (HAT) segment separate from the pelvic segment. Additionally, the upper part of the body may also be modelled to include segments for the lumbar region, thorax, arm, forearm, hand, neck and head segments. Body segments are defined in terms of their physical dimensions, mass, location of centre of mass and moments of inertia. Each body segment also includes its own Cartesian coordinate system, usually aligned with standard anatomical directions of the segment (e.g. proximodistal, anteroposterior, mediolateral). Other coordinate systems can be defined relative to that of the segment or to the global system and allow for ease in defining model geometry. The segments further rely on the assumption that they do not deform during locomotion (i.e. are rigid). In the software modelling/simulation environment, body segments are almost universally represented by virtual surface models of the relevant skeletal element(s).
Figure 1. MoCap model mid-simulation of walking. The image shows the model (rigid body segments represented by skeletal elements and visual representation of muscle models), four force plates (grey), the ground reaction force vector (blue line through model), force plate readings (blue lines below the force plate), experimental marker locations (blue spheres), virtual markers (red spheres with coordinate system arrows) and global coordinate system (yellow arrows).
Body segments are connected to other segments at joints, which represent articulations (e.g. hip, knee, ankle) within the body. Joints are frictionless and non-deformable and, consequently, forces and moments transfer across these joints in an equal and opposite manner. A joint definition includes its location relative to the body segments it connects, as well as the degrees of freedom (DOFs) permissible at the joint (one, two or three axes of rotations and, less commonly, translations). The direction of rotation axes are usually expressed in terms of a local (body segment) coordinate system.
With the body segments and joints defined, the skeletal portion of the musculoskeletal model is sufficient to conduct dynamic analyses. Musculoskeletal models extend the analysis by including models of individual muscles or even portions of muscles (figure 1). Muscle models range in their complexity and, as with other parts of the model, the necessary complexity depends on the motion being simulated. (For a description of how muscle complexity can be modelled, see §5.3.) A muscle model minimally requires geometry (i.e. path definition) and a parameter representing its maximum possible force. More complex muscle models also include parameters for pennation angle, optimal fibre length, physiological cross-sectional area (PCSA) and activation and deactivation times. Muscle geometry must include the coordinates of its origin (location relative to the proximal segment) and insertion (location relative to the distal segment), but muscles can be geometrically more complex than a straight line path between origin and insertion. The most basic of these complexities derives from a muscle with a broad origin (e.g. gluteus medius) or insertion (e.g. adductor magnus). In these muscles, a single origin-to-insertion path does not encompass their possible actions. One way to solve this is to split the muscle into muscle parts, referred to here as muscle elements. Each muscle element represents a region of the muscle (e.g. the superior portion of adductor magnus) and can have its own set of muscle parameters.
Most modelling systems allow for additional coordinate positions that are used to define the curved path that the physical muscles take between bony attachments (e.g. the path sartorius takes across the anterior thigh as it traverses from anterior superior iliac spine (ASIS) to proximal tibia). Via points, the simplest type of path point, are one or more coordinate positions that are intermediate between the origins and insertion points. The muscle assumes a series of straight line paths between the origin, via point(s) and insertion. The location of the via point is fixed relative to the segment's centre of mass and the connection to other segments. With via points, the muscle is required to pass through each specific point at all times. These points allow muscles to pass over other body surfaces (e.g. iliacus wrapping over the pelvic brim). In addition, muscle geometry paths can be altered by the use of ‘wrapping surfaces’ that are designed to account for relative motion between musculoskeletal elements and other soft tissue (e.g. rectus femoris sliding across the soft tissue structures anterior to the hip). These wrapping surfaces are geometric objects (e.g. cylinders, ellipsoids) that are fixed relative to a particular body segment. Muscle paths are not constrained to pass through any particular coordinate point (as with via points), but rather can move freely over the wrapping surface to find the shortest possible path across the surface between origin, insertion or via points that are beyond the wrapping surface. As the proximal and distal segments change orientation relative to each other, the point of contact between the muscle and the adjacent musculoskeletal elements can move (e.g. gastrocnemius relative to the posterodistal femur as the knee flexes). Forces developed from the action of the muscle act on the body segment used to define the geometry (origin, insertion, via points and/or points on wrapping surfaces). Muscles may also be modelled to be more complex in terms of the number of parameters that define force generation. Minimally, a muscle model must include a maximum force value which defines the force when the muscle is maximally activated, but other approximations of muscle force generation are available, e.g. the Hill-type muscle model [64].
The final required model element is a set of virtual motion markers, which is needed to match model motion to motion measured experimentally (figure 1). For motion tracking (described below), the human subject is fitted with a set of markers that are tracked by a camera system. The musculoskeletal model must include one virtual marker for each experimentally measured marker. Each virtual marker is connected to a body segment and defined in terms of that segment's coordinate system. For instance, if researchers affix a motion-tracking marker over the ASIS of a human participant, the model will require a virtual marker connected to the pelvis body segment with coordinates that position it correctly relative to the model ASIS.
2.2. The input data: motion data and ground reaction forces
An inverse dynamics approach to musculoskeletal modelling (described below) requires that the motions (translation and rotation) of the rigid segments (i.e. the skeletal elements) and externally applied forces be known. This process is well established and has been described in detail elsewhere (e.g. [65–67]). To obtain motions, retroreflective (infrared) markers attached to the skin or clothing are frequently used. Each body segment requires a minimum of three markers that are trackable throughout the entire sequence of movements, so this traceability should be checked in all portions of the data collection volume. The three-dimensional ground reaction force (GRF) for a complete stance is also necessary. GRFs are typically acquired through force plates embedded either in the floor or in a treadmill.
2.3. The simulation
2.3.1. Scaling and marker optimization
Model scaling and virtual marker position optimization are the processes of determining the model segment dimensions and virtual marker positions so that the virtual model matches the human subject and experimental protocol (motion marker placement on the human subject) as well as possible [68]. The degree of matching between the model and experimental data is evaluated using a global error metric between experimental and virtual motion marker positions. During the inverse kinematic analysis (described below), musculoskeletal modelling software finds the position of the model, including overall location in space and individual joint positions, that minimizes the sum of squared deviations between the virtual marker set and experimental motion marker data. The goal of marker optimization and body scaling is to minimize this error term.
The human models that are included with many modelling software packages are generic human models. For analysis, the model must match the subject for whom the input data are known in terms of mass, limb parameters (e.g. distances between joint centres and segment moments of inertia) and virtual marker positions [69,70]. Marker position optimization and scaling are intertwined processes because marker positions can, in part, determine segment dimensions. Initial scaling of the generic virtual model achieves a loose match between the overall bodily dimensions of the subject and those of the virtual model. If a virtual model from a repository is used (rather than creating a virtual model from scratch), then the model will have been developed using the parameters of the individual whose data were available (i.e. a generic model). This generic model will almost certainly not represent the dimensions of the subject on whom motion-tracking data and GRFs were collected, so the generic model needs to be scaled [69]. The simplest scaling approach is uniform scaling of lengths and volumes. Total body mass is the most basic dimension needed, and the ratio of the body mass of the subject to that of the generic model is typically used to adjust segment masses and other parameters that are volume-dependent. The ratio of subject to generic model stature, a length scale, is used as an initial estimation of segment dimensions and other dimensions that are length-dependent. Other anthropometric data (e.g. distances between landmarks) can be used instead of, or in addition to, stature. Other scaling methods are also available, including fully manual scaling, using body fat to adjust mass distribution among the segments, and hybrid scaling (e.g. manually scaling some parameters and algorithmically scaling others). Furthermore, depending on the sophistication of the software used to complete the mechanical analysis, these initial scalings may be adjusted during marker optimization.
Different software packages conduct marker optimization and secondary scaling of segments using different algorithms and require varying levels of interaction/input from the user. At one end of the user-input spectrum, some systems require that the user optimizes the locations of the markers on the model ‘by hand’, moving them relative to their parent segment and then interrogating error values associated with each marker during motion tracking. This process is tedious and may require significant effort before optimal marker position is achieved. At the other end of the spectrum, some software packages change both segment dimensions and virtual motion-tracking markers automatically in a single cohesive optimization process that requires little user interaction. In these frameworks, the user can permit or restrict changes in virtual marker position in none, some or all of the coordinate directions. For instance, motion-tracking markers associated with easily palpable bony landmarks (e.g. medial malleolus) are also easy to locate on the virtual model (referencing the virtual skeleton) and may be considered ‘fixed’ in their location relative to their body segment (the shank). Thus, during the optimization process, such a fixed marker would not change its position relative to the segment, but would influence the determination of shank length. Other markers (e.g. marker plates attached mid-segment) may be allowed to optimize their position in all three coordinate directions, and then contribute little to segment length determination. Finally, some markers may be well defined in two coordinate directions, but less so in the third. For instance, the superoinferior and mediolateral positions of an ASIS marker may be well defined and constrained from moving relative to the pelvis during the optimization process. The anteroposterior position of the virtual ASIS marker relative to the virtual pelvis may be more difficult to estimate as it reflects tissue thickness over this bony landmark (which is hard to measure on living subjects in the absence of internal imaging modalities). Thus, the user might fix the ASIS marker in two coordinate directions, but allow its position to optimize in the third. The fixed mediolateral and superoinferior position on the two ASISs would then influence the scaling of pelvic width.
2.3.2. Inverse kinematics/motion tracking
To carry out inverse dynamics and, in turn, estimate internal forces and moments, segment linear and rotational accelerations must be obtained. Most musculoskeletal modelling software (e.g. Anybody Modeling System™, OpenSim) uses a kinematic analysis procedure of over-determined biomechanical systems, which is formulated as a weighted least-squares optimization problem that matches experimental and model marker positions (e.g. [47,71]).
2.3.3. Inverse dynamics
While statics is the branch of mechanics that deals with stationary objects, dynamics is the branch in which motions are the focus. In statics, an equilibrium that greatly facilitates examination of the system exists where the environment (e.g. the substrate) and the structure (e.g. the human body) are balanced. In standing, the body pushes through the feet and against the ground exactly as much as the ground pushes against the body through the feet. In dynamics, the balance is time-dependent and the solution to the system is through the equations of motion, which relate forces to motions. When the forces that produce the motions are fully known, the problem is characterized as forward dynamics, while when motions are known it is described as inverse dynamics.
In human bodies, the measurement of motions (via marker tracking or other means) is much less invasive and more feasible than measuring the muscle forces that produce the motions, so currently most analyses that depend on the musculoskeletal forces generated in human gait are typically treated as inverse dynamic problems (e.g. [72–74]). Further complicating the solution in multi-body problems (i.e. those with multiple, linked rigid bodies) is that the system can be interrogated from an external or internal perspective. The external view includes the action of the GRF to produce the movement of the body along a path, while the internal one includes both the joint reactions and the action of muscle forces in creating changes in joint angles. In mechanics, the relationship between the number of unknown variables and the number of known variables determines the difficulty of the solution. Fortunately, the internal and external perspectives can be exploited to develop solutions. For instance, in the single-stance foot, the GRF is the external force and produces the moment acting on the foot segment and reacted to by the ankle joint reaction force and moment (which are an internal force and moment). Once the ankle joint reaction force and moment are known, they can be combined with the motion of the shank to solve for the knee joint reaction and moment. The inverse dynamic solution allows for this time-dependent chain solution to determine the joint reactions.
2.3.4. Muscle redundancy solution
A system where the number of variables with unknown value equals the number of independent equations that fully define the system is determinant and a unique solution exists. For systems where there are more variables with unknown values than equations—called indeterminate—assumptions are made to provide more equations [75,76]. Many biological problems, including developing muscle forces, are indeterminate because there are more muscles capable of acting across a joint than necessary from a purely mechanical perspective, i.e. some muscles are (mechanically) redundant. If the desired result from an analysis that uses inverse dynamics is the development of muscle forces, then the joint moments determined from the inverse dynamic solution need to be apportioned to individual muscle forces. This is not a trivial undertaking, because the degree of redundancy is often quite large, so no general consensus on the most accurate criterion on which to base an apportionment algorithm currently exists [77].
Musculoskeletal modelling solutions typically include an algorithm that uses muscle parameters to apportion the joint moments to particular muscles (as muscle forces). Early solutions used PCSA to solve the muscle redundancy problem [75]: the apportionment was made based on the ratio of the PCSA of a particular muscle relative to the total PCSA of all muscles capable of acting (e.g. [61]). While effective from an algorithmic perspective, PCSA-based apportionment does not reflect the natural solution, because cross-sectional area approaches assume equal activation of all possible muscles.
Another algorithm for solving the muscle redundancy problem attempts to minimize the sum of muscle activations raised to a power [76]. The activation represents a proportion of total maximum strength the muscle can generate given the current state of the model (i.e. activation is generally bound between 0 and 1). If the muscle activation exponent is 1, i.e. a linear combination, force is allocated to the muscle with the greatest moment arm until it reaches its maximum allowed value, and then other muscles are recruited. Exponents greater than 1 will distribute forces more evenly across muscles available to produce the moment at a particular joint. The greater the power, the more even will be the distribution (PCSA-based activations are essentially an exponent of infinity). Other criteria for apportionment of the joint moments have also been advanced, such as electromyography (e.g. [78,79]), but activation-based methods remain the standard [76]. Until in vivo muscle forces can be reliably measured, all solutions to the muscle redundancy problem should be treated as hypotheses about how joints and muscles produce motion.
2.4. The output
A valuable aspect of musculoskeletal modelling (or of any model) is that variables of interest can be interrogated anywhere within the model. Standard output from a musculoskeletal model includes linear and angular displacements, velocities and accelerations for all body segments. Similarly, joint positions, velocities and accelerations are also tracked and can be exported for analysis. Muscle length, activation, moment arm and force are all possible outputs, as are joint reaction forces. Additionally, users can track the position of muscle origins/via points/insertions throughout the modelled motion. Conveniently, joint reaction forces can be output in the global coordinate system or the coordinate system of the relevant body segments.
Human locomotion requires the repetition of cycles of limb and body motions that can be reduced to a stride. For walking, a step is defined from an event (e.g. initial contact) on one foot to that same event on the contralateral foot and a stride is two contiguous steps. A stride can vary in duration, even at a constant velocity, so creating an average stride or evaluating the variability of stride parameters requires normalization of the strides to some consistent duration. The typical choice is per cent of stride cycle, where the stride duration is described in terms of equal portions of the entire stride (e.g. [80,81]). Evaluating the variability of parameters of interest (e.g. muscle force, joint reactions) can be an important step in validating the calibration of a solution. While some differences among trials of the same individual are expected, an evaluation of the variability in light of the question that the model is designed to answer should occur.
2.5. Human variability
A crucial advantage of musculoskeletal modelling is the opportunity to include aspects of morphological variation, including the production of subject-specific models. Model scaling and motion tracking incorporate some aspects of human variability, including differences in body mass, limb lengths/dimensions and kinematics, but other aspects are less obviously accessible for modification in musculoskeletal modelling software. These less accessible aspects, such as variation in bone or muscle morphology, are also important to consider. Quantifying human morphological (musculoskeletal) variation is at the heart of biological anthropology. For example, human (and hominin) pelvic morphology has a profound influence on bipedal capabilities [14,80], and has been the focus of intense scrutiny [41,81,82]. Humans are known to be sexually dimorphic in pelvic size and shape [83–85]. Furthermore, human pelvic morphology is also known to vary by population [86–88]. Sources of variation can be incorporated into musculoskeletal modelling to understand pathology [31] and differences in joint reaction forces [35]. For example, marker-informed parameter optimization and secondary scaling (described above) can change the width of the pelvis if experimental and virtual markers are included on bony landmarks of the pelvis, such as the ASISs. If the subject is wider than the generic model, the parameterization can widen the pelvis so that the virtual and experimental markers are more closely aligned. Scaling can also be accomplished manually either by applying factors to the three segment coordinate directions or by redefining the locations of specific anatomical landmarks associated with the segment. The former method moves all points proportionately while the latter allows for localized warping.
In addition to skeletal variation, aspects of muscle morphology that vary among humans may also be critical for accurate musculoskeletal modelling. For instance, variability in muscle parameters has been of concern for some time [89]. While repository models represent a particular person's morphology, people differ in important musculoskeletal modelling parameters such as muscle volume, optimal fibre length and pennation angle (table 1).
Table 1. Comparison of muscle parameters between the standard MoCap model in AMMR™ v. 2.3 (AnyBody Technology, Aalborg, Denmark; from Klein Horsman et al. [90] and Carbone et al. [34]) and Charles et al. [91]. Muscle names are per the MoCap AMMR v. 2.3. Where multiple muscle regions are present in the MoCap model, the values for the muscle volumes are summed while the optimal fibre lengths and pennation angles are averaged.Collapsemuscleoptimal fibre length (cm)pennation angle (°)muscle volume (ml)MoCapCharles—averageMoCapCharles—averageMoCapCharles—average
adductor brevis | 10.38 | 7.6 | 0.0 | 11 | 108.90 | 93 |
adductor longus | 10.57 | 11.0 | 0.0 | 12 | 159.56 | 159 |
adductor magnus | 9.97 | 23.1 | 0.0 | 12 | 559.25 | 567 |
biceps femoris caput breve | 9.14 | 10.9 | 0.0 | 9 | 107.95 | 92 |
biceps femoris caput longum | 8.54 | 20.4 | 29.9 | 11 | 232.01 | 194 |
extensor digitorum longus | 5.99 | 13.8 | 8.3 | 7 | 32.29 | 76 |
extensor hallucis longus | 5.98 | 10.6 | 14.4 | 7 | 36.27 | 21 |
flexor digitorum longus | 3.84 | 28.4 | 25.28 | |||
flexor hallucis longus | 5.20 | 30.1 | 161.50 | |||
gastrocnemius lateralis | 5.69 | 12.2 | 25.4 | 9 | 136.36 | 128 |
gastrocnemius medialis | 6.01 | 9.7 | 10.8 | 10 | 263.26 | 230 |
gemellus | 3.43 | 0.0 | 28.40 | |||
gluteus maximus | 13.57 | 0.0 | 936.55 | |||
gluteus medius | 6.76 | 5.3 | 674.71 | |||
gluteus minimus | 3.28 | 0.0 | 82.59 | |||
gracilis | 18.11 | 17.3 | 0.0 | 6 | 87.97 | 91 |
iliacus | 8.16 | 8.8 | 203.13 | |||
obturator externus inferior | 6.88 | 0.0 | 37.88 | |||
obturator externus superior | 2.77 | 0.0 | 68.18 | |||
obturator internus | 2.05 | 0.0 | 52.08 | |||
pectineus | 9.50 | 0.0 | 64.58 | |||
peroneus brevis | 4.54 | 23.1 | 86.09 | |||
peroneus longus | 5.08 | 15.8 | 121.52 | |||
peroneus tertius | 4.27 | 19.1 | 26.52 | |||
piriformis | 3.88 | 0.0 | 31.25 | |||
plantaris | 4.77 | 0.0 | 11.36 | |||
popliteus | 2.40 | 7.4 | 0.0 | 8 | 25.57 | 15 |
psoas major | 9.92 | 13.4 | 193.18 | |||
psoas minor | 5.91 | 0.0 | 6.63 | |||
quadratus femoris | 3.37 | 0.0 | 49.24 | |||
rectus femoris | 7.84 | 14.2 | 21.9 | 8 | 226.33 | 249 |
sartorius | 34.71 | 40.8 | 0.0 | 205.49 | 143 | |
semimembranosus | 8.09 | 15.8 | 25.0 | 12 | 138.26 | 244 |
semitendinosus | 14.16 | 18.3 | 0.0 | 8 | 208.33 | 186 |
soleus | 4.40 | 14.6 | 61.6 | 12 | 792.91 | 461 |
tensor fasciae latae | 9.49 | 0.0 | 83.33 | |||
tibialis anterior | 6.83 | 13.7 | 9.6 | 7 | 181.49 | 129 |
tibialis posterior | 3.78 | 34.2 | 163.36 | |||
vastus intermedius | 7.68 | 18.1 | 11.8 | 11 | 292.61 | 521 |
vastus lateralis | 6.69 | 19.6 | 0.0 | 15 | 583.33 | 606 |
vastus medialis | 7.16 | 15.9 | 0.0 | 14 | 454.27 | 415 |
Several possibilities for future work that explores the impact of morphological variation on muscle and joint forces are worth noting. For instance, while some anthropological work has elucidated the interactions of morphology and walking conditions (e.g. velocity, changes of direction, burdens and gradients) on energy expenditure (e.g. [92–94]), much less is known about how morphology and walking condition factors interact to affect joint and muscle forces (although see [31] for an example). People vary morphologically and walking under varied conditions is an activity of daily living, so this area is an important one for many disciplines to pursue including biological anthropology, biomedical engineering and product design. Perhaps more importantly, collecting the requisite data to inform a musculoskeletal model is time consuming, invasive and dependent on (at this time) fragile laboratory equipment. Consequently, acquiring samples that represent the breadth of human morphological variation is unlikely to occur in the near future. Once validated, though, a musculoskeletal model can be modified to assess the impact of these known morphological differences on muscle and joint forces. Extending that line of research further, a longer-term goal would be use musculoskeletal modelling to understand the influence of the morphological differences between modern humans and extinct hominins on muscle and joint forces, an area that has sought biomechanical solutions without effective software support for decades [12,14,61,95,96].
2.6. Current study
Here, we demonstrate the process and output of one human subject during walking using AnyBody Technology's musculoskeletal modelling software (Anybody Modeling System™; AnyBody Technology, Aalborg, Denmark). We present several of the possible outputs that could be of interest to biological anthropologists. This includes joint kinematics, joint reaction forces and a selection of lower limb muscle force magnitudes. We make some recommendations for future researchers to help facilitate their use of software and output. Additionally, by presenting a set of possible outputs, we demonstrate the motivation for undertaking musculoskeletal modelling and that results may have a direct bearing on anthropological questions, or as intermediate steps to other analyses (e.g. finite-element analysis). Additionally, while engineering has much to offer biological anthropologists, the latter have an appreciation for human variation which may enhance future generations of musculoskeletal models.
3. Material and methods
3.1. Participant
One healthy participant was used for this demonstration, but he is part of a larger study aimed at examining human walking. The participant was 36 years old (male; body mass, 74.3 kg; stature, 172 cm) and reported being free from lower-limb injuries. The University of Washington's Institutional Review Board approved all aspects of this study (IRB no. STUDY00001125).
3.2. Experimental protocol
Kinetic and kinematic data were measured using a four-force plate (Kistler, Switzerland), 10-camera motion capture system (Qualisys, Sweden) in the Amplifying Movement & Performance (AMP) laboratory of the University of Washington. Thirty infrared-retroreflective 5 mm hemispherical markers were affixed to anatomical landmarks and used to track motion through the laboratory (figure 1; electronic supplementary material). The participant walked unshod the length (10 m) of the gait laboratory five times at his self-selected preferred walking speed (1.15 m s−1). The participant walked in a straight path with his eyes directed at a point on the far wall of the laboratory and contacted the surface of each force plate with a single foot. Trials with multiple feet on a single force plate or ones where the stance foot crossed plates were immediately redone. The participant was allowed several attempts prior to data collection to become familiar with the protocol. All trials exhibited similar kinematics and kinetics (data not shown), so we selected one trial to use in this demonstration.
Marker data and GRFs were collected at 120 Hz and 1200 Hz, respectively, and then were filtered at 5 Hz and 10 Hz, respectively, with a fourth-order, low-pass zero-lag Butterworth filter. Calibration of the system yielded a limitation in its fidelity for marker data of 1 mm and force data of ±2.5 N for the direction of travel (X), ±5 N side (Y) and ±25 N vertical (Z). All data were exported from the Qualisys Track Manager software in C3D file format, which can be read directly by the modelling software (provided in electronic supplementary material).
3.3. Musculoskeletal model
We used the ‘MoCap’ model from the AnyBody Managed Model Repository™ (AMMR v. 2.3; AnyBody Technology, Aalborg, Denmark; http://www.anybodytech.com/software/model-repository-ammr), which is a multi-trial, full-body, motion-capture-driven human gait model, and the commercially available AnyBody Modelling System™ (v. 7.3; AnyBody Technology, Aalborg, Denmark) to calculate forces in the lower limb [34,97]. The MoCap model that we used for this demonstration has been used to assess human gait in numerous applications (e.g. [98,99]). The MoCap model is assembled from a trunk and left and right lower limb components. The trunk includes a lumbar region, trunk, neck and head. The lower limbs consist of the following segments: thigh, patella, shank, talus and foot. Each lower limb has six total DOFs, including all three rotations at the hip and one each at the knee (flexion/extension), ankle (plantar flexion/dorsiflexion) and subtalar (inversion/eversion) joint. The pelvis (relative to the ground) has six DOFs, three translational and three rotational. Forty-one anatomically defined, bilateral muscles (table 1) are represented by 169 muscle elements in each model lower limb (e.g. gluteus medius is composed of 12 separate muscle element actuators) [34]. The muscle elements in this model generate force as a function of their activation and a strength parameter which is appropriate for human walking [97,100]. We defined a virtual motion-tracking marker for each of the experimentally measured motion markers (electronic supplementary material). The polynomial algorithm that was employed to solve the muscle redundancy problem apportioned force to muscles using an objective function with power = 3. (For more information regarding the model, see §5.3.)
3.4. Musculoskeletal simulation
The first step of the simulation was to scale the model to match the dimensions of the subject while simultaneously optimizing the marker coordinate and joint rotation axes' locations [68]. After initial scaling, the values of these parameters (i.e. segment dimensions, marker locations, rotation axes' locations) were determined using the marker motion data from a trial. Values for all parameters were optimized to minimize the sum of square distances between model marker positions and experimental marker positions. Following convergence of parameter optimization, we conducted an inverse kinematics analysis. Finally, we conducted the inverse dynamic analysis, which includes apportioning muscle forces using the algorithm that minimizes the cost function as a sum of muscle activation values raised to the third power. This value has been shown to allocate muscle activations in a way that reflects experimental data [77,100,101].
3.5. Output variables and data processing
We output several kinematic and kinetic variables, including joint kinematics, filtered GRF components and joint reaction forces (electronic supplementary material). Additionally, we output the force magnitude for a selection of the muscle elements for the left lower limb during an entire stride cycle of walking. Muscle element forces within anatomically defined muscles were summed for presentation (e.g. 12 elements that represent gluteus medius). All time-related curves were resampled to 1% increments of the stride cycle.
4. Results
GRF component curves across the gait cycle are presented in figure 2. As is typical for human walking, the vertical component has two peaks (early stance peak, 719 N; late stance peak, 778 N) on either side of a valley (633 N) which occurs around the middle of stance phase. The average of the two peaks is 103% of body weight (729 N), while the valley represents 87% of body weight.
Figure 2. Filtered ground reaction force components. Grey vertical line represents moment of toe-off.
Pelvic rotation and obliquity and hip, knee and ankle kinematics are provided in figure 3. All kinematics are of the left lower limb throughout a single stride. Pelvic motion is described in terms of the right hip (previous stance limb) relative to the left hip (current stance limb) [102]. Hip, knee and ankle joint reaction forces, which include contributions from muscle forces, are provided in figure 4. In each joint, the vertical joint reaction force is greater in magnitude than the anteroposterior and mediolateral components. All three joints show a large peak during the second half of the stance phase, and this late peak is significantly larger in the ankle. Muscle force profiles for 18 anatomically defined muscles across the stride cycle are provided in figure 5. During the stance phase of walking, gastrocnemius, soleus and gluteus medius generate the largest forces. All output data are provided in the electronic supplementary material.
Figure 3. Angular kinematics for the pelvic segment and hip, knee and ankle joints. Vertical grey line represents toe-off.
Figure 4. Joint reaction forces expressed in the global coordinate system. Vertical grey line represents toe-off.
Figure 5. Muscle forces for 18 anatomically defined muscles.
5. Discussion
5.1. Model output
The subject that we used to demonstrate the musculoskeletal modelling approach to modelling walking is an adult male without gait pathology and within species-typical anthropometric characteristics. Consequently, the model-derived motions should accord with those determined experimentally and extensively described in the last 40 years of gait research (e.g. [66,103–105]). The model results presented above do meet this requirement, as discussed below in detail. We focus on the model-derived kinematics, because these kinematics are the foundation for all subsequent analyses in musculoskeletal models. Any analysis should begin by verifying that model-derived variables align with those obtained empirically, and the cause of any discrepancies should be understood or noted as a potential limitation. We belabour the details in the following paragraphs in order to demonstrate the basic approach to model validation.
In our subject at heel strike, the pelvis is rotated approximately 10° in the horizontal plane such that the right hip is posterior to the left hip (e.g. the current stance limb). The pelvis then swings through a neutral position and, while leading up to the next right heel strike, into an orientation with the right hip anterior to left hip (approx. 10°). The pelvis drops in the frontal plane with the right hip inferior to the left (stance) hip soon after heel strike as the right lower limb moves into the swing phase with pelvic obliquity ranging between ± 4°. The pelvis remains in this position until the subsequent right heel strike, after which the pelvis approaches the neutral position. Model-derived pelvic rotation and obliquity are similar to empirically determined values (e.g. [103]).
In the sagittal plane, the left hip begins the stance phase in flexion of approximately 17° and moves through a neutral position and into extension of approximately 17°. Winter [66] describes a less symmetrical flexion–extension cycle at natural cadence (mean values of approx. 20° of flexion and 11° of extension), but the model-derived hip rotations are well within one standard deviation of his mean values. During most of the stance phase, the left thigh is adducted less than 6° in accord with others (e.g. [103]).
While the pelvic and hip joints have three rotational DOFs, the knee and the ankle only rotate in the sagittal plane. The left knee flexes slightly (10°) early in the stance phase, but maintains a near fully extended position through midstance, after which it begins to flex in preparation for the swing phase where flexion exceeds 60°. After the initial heel strike and as the forefoot contacts the substrate, the ankle plantar flexes slightly (less than 5°) and then assumes a dorsiflexed position for most of stance. Rapid plantar flexion occurs in late stance as the body is propelled forwards, reaching a plantar flexion value of greater than 10°. The ankle returns to a dorsiflexed position for swing. These values for both knee and ankle angles are typical (e.g. [66,103,105]).
Consequently, the kinematics of the model match those that we expect empirically. While this check is not sufficient to ensure that the internal (muscle and joint reaction forces) are reasonable, it is a necessary condition and validation for any model. Incorrect kinematics guarantee incorrect or suspect muscle forces.
Joint reaction forces have been empirically measured only with instrumented joint implants. By definition, the surgery creates conditions that limit applicability of the results, but no other method currently exists to measure force. Pedersen et al. [45] measured hip joint contact forces in a 72-year-old man (63.2 kg) with a hip replacement and found that the vertical force was 1750 N, considerably less than the almost 3500 N we found with the musculoskeletal model. Important differences between this patient and our subject exist, including age (72 versus 36 years), body mass (62.3 versus 74.3 kg), health status (58 days post hip replacement surgery versus healthy) and walking speed (0.89 versus 1.15 m s−1). While hip and knee joint reaction forces presented here are larger than in vivo values using instrumented prostheses [77,106,107], our model-derived joint forces are, however, similar to other models of healthy adults [35,107].
Although it is difficult to compare muscle forces with experimental data, it is possible to consider the pattern of activation with reference to the muscle's action. The muscles that generate the highest force in the model are gluteus medius, soleus and gastrocnemius. Gluteus medius plays a critical role in bipedal walking by preventing pelvic drop during single stance [103,105]. The gluteus medius muscle profile from the model is double peaked, and the timing of the two peaks is coincident with the two peaks in the vertical GRF, indicating the muscle's support of the pelvis during single stance. Gastrocnemius and soleus plantar flex the ankle; when the foot is in contact with the ground and the hip is extended, plantar flexing the ankle moves the body anteriorly [103,105]. The gastrocnemius and soleus muscles show the greatest force generation during the propulsive portion (second half) of the stance phase, in keeping with our expectation from their function.
Other muscles also display moderately high forces. Tibialis anterior, which dorsiflexes the ankle, shows muscle force generation immediately after heel strike when the body begins to move over the foot, and then is mostly inactive until the second half of the swing phase. The quadriceps demonstrate a force peak early in stance phase, and again at the stance–swing transition, potentially relating to their role in extending and stabilizing the knee [105]. The hamstring muscles show the greatest peak at the end of swing phase and into early stance phase. Activation of these muscles at the end of swing phase is associated with decelerating the limb prior to heel strike. Early in stance the GRF is directed anterior to the hip and knee, and the hamstring muscles are active to prevent further flexion (hip) and extension (knee) of these joints. The hip flexors are active late in the stance phase, and may be active to counteract hamstring hip extension moment as the hamstrings flex the knee in preparation for swing phase. The adductors show initial peak in stance (adductor magnus) and late stance (adductor longus), but generally smaller than forces in other muscles.
5.2. Input data choices
The model-derived motions and joint reaction and muscle forces align with our expectations from the experimental data in the literature. We achieved these results through careful attention to detail when we collected the data we used as inputs to the musculoskeletal model. We provide below some brief suggestions from our experience. Because motions are derived from the change in position of these markers, marker placement on the segment should be chosen to define fully the motion of the segment and to reduce marker positional noise as much as possible. Positional noise arises from the inherent error of measurement of the marker's position in space (system error) and from movement of the surface (skin or clothing), where the markers are placed relative to the underlying rigid segment (e.g. [108–110]). System error is reduced through calibration; for our system, it is less than 1 mm for any marker. Movement between the external (skin) surface and the bone is more complicated to reduce, but may not be clinically relevant [110].
Marker position is a compromise with soft tissue coverage of underlying bone, landmark palpability (to aid placement) and distance between markers (such that marker positional error is much less than the distance between markers) all being important. The shank and foot have palpable landmarks that are covered by minimal soft tissue (e.g. the malleoli), but the thigh (e.g. greater trochanter) and, especially, the pelvis are more challenging segments on which to affix markers. This problem is exacerbated when the subjects of interest have more soft tissue (muscle or fat) overlying the bony landmark, and considerable soft tissue coverage is common in the hip–pelvis region. One approach to ameliorate some of this difficulty is to use a marker plate which has four (or more) markers mounted on a hard plastic plate that is then strapped to the segment. If a plate can be affixed firmly to the segment, as it can be to the thigh, the marker plate eliminates between-marker error. Unfortunately, marker plates are difficult to affix firmly to the pelvis. Marker placement issues, consequently, are most difficult to overcome to track the motion of the pelvis. Yet, for musculoskeletal models of human locomotion, pelvic motions are critical components of system performance. We achieved the motions shown in figure 3 through extensive protocol development and testing.
Another critical input to the simulation are the GRFs. It is possible to conduct an inverse dynamics solution of the single support phase of walking with a single force plate (e.g. [111]), but extending the analysis to the double support portions of the stride is more difficult with one force plate because it would require estimation of the unknown contribution of the other stance foot. Consequently, any research that requires information about double support should optimally be conducted with multiple force plates. The placement of the force plates relative to each other should be appropriate to the population and gait. Stride lengths for adults or running gaits are longer than for children or walking, and plates should be positioned to prevent subjects ‘targeting’ foot placement or large numbers of repeated trials. Data analysis is less complicated if stance occurs on one plate, but the data from two immediately adjacent plates can be used as long as there is no gap between them [112]. Multiple feet contacting a single plate is also possible to resolve with additional information (such as localized pressure sensors) or predictive methods (e.g. [113,114]). All of these remediation techniques are time-intensive. The laboratory in which the presented data were collected has four floor-embedded force plates that are positioned to allow normal stride lengths for adults when walking. The subject and trial that we used for this demonstration is part of a much larger study [115]. Nonetheless, we assessed each foot placement on all four plates for each trial before moving on to the next trial. Such data checking is a tedious effort at the time of data collection, but well worth the effort to ensure all data are suitable for future analyses. Using four force plates, we had data frames before and after the data that we used in the analysis. This allowed us to discard the initiation and termination frames of the simulation, which are frequently unusable. Consequently, researchers are advised to locate force plates well and test the placement prior to data collection.
Finally, researchers may consider collecting complete anthropometrics of subjects and using these in the initial scaling as the starting points for the parameterization. People vary considerably in terms of their proportions [35,37], so standard models, such as the model we used here, represent a morphology that may not match the dimensions of any particular subject. If your model starts with segment lengths that are substantially different from the segment lengths of the model, the parameterization algorithm may not be able to converge on a suitable solution. Two failures can occur: the parameterization can fail to converge, which produces an obvious error message, or the parameterization converges, but, in order to converge, the simulation requires changes to the morphology or motions that are inconsistent with the experimental data or subject dimensions. This latter failure can be insidious. An absolutely critical step in model validation is to check the motions and shapes of the segments after parameterization. A poorly parameterized model is incorrect and will not produce acceptable motion or joint and muscle forces.
5.3. Modelling choices
In order to produce this demonstration, we made a number of choices with regard to the model which bear discussion. The most important of these was the choice to use a standard, repository model, in this case the MoCap model of AMMR v. 2.3. By using the MoCap model we accepted many assumptions and compromises, but gained the experience and efforts of many dedicated musculoskeletal modellers. Nonetheless, the responsibility of assessing the ability of any model to answer the question of interest rests with the researchers using the model. Our intention was to demonstrate some of the decisions and validations needed in musculoskeletal modelling of human walking and our choices reflect that.
The first of these choices is the theory on which the muscle actuation rests. The muscle models used here generate force as a function of maximum force and activation. More sophisticated muscle models, such as the Hill muscle model, are available. The Hill muscle models are thought to provide the most complete model of muscle behaviour [116], but the Hill model requires several parameters for its full definition, including: maximum isometric force, tendon slack length, optimal fibre length, pennation angle at optimal fibre length, activation and deactivation times, as well as maximum contraction velocity. The number of parameters that must be defined for a Hill muscle model, as well as the sensitivity of muscle behaviour to parameter values, means that considerable effort must be expended to parameterize and optimize Hill-type models. As with all models, poorly informed inputs (i.e. muscle parameters) lead to unstable results. This means that for certain activities it may be advisable to use a simpler, robust model rather than a more complex, unreliable one. For human activities during which muscles operate at slow contraction velocities and near their optimal fibre length (e.g. human walking), simple muscle models that depend only on geometry and a maximum force parameter are preferable [100]. Consequently, we accepted the use of the simple muscle model employed in the MoCap gait model.
Another choice was the selection of the scaling method. Understanding which scaling method to use (e.g. mass–length–fat), requires understanding the underlying assumptions of the method (e.g. distribution of mass to segments), to which model parameters the scaling applies (e.g. which segment length scale to which muscle length) and how the model uses the scaling to define parameters (e.g. are all parameters scaled or just a subset?). We made many other choices by not changing components that are inherent to the generic MoCap model, including acceptance of the form of the algorithm used to solve the muscle redundancy problem. Such choices are often embedded in standardized models, and the user does not have to modify any of them in order for the model to run to completion. This makes ensuring that the question of interest can be addressed by the model become an intentional assessment. For instance, the standard MoCap model has only one rotational DOF for the knee, which works well to assess overall walking parameters but is less well suited to a detailed study of the knee. Furthermore, standardized models are improved as more experience and data become available (e.g. [97]). For instance, the MoCap model's leg (TLEM v. 2.1) has been modified to reflect muscle paths more accurately by the inclusion of wrapping surfaces.
We also relied on the standardized definition of available outcome variables such as muscle forces and joint motions. Given that most software allows for forces and motions to be expressed in different ways (at a minimum, in different coordinate systems) it is imperative to understand the context of the reported variable. Is an output reported in the global or segment coordinate system? Is a joint force an internal force or a reaction? In what direction does the muscle force act? As with all big data, extensive data exploration and a complete understanding of the theories (mechanics), procedures (inverse dynamics) and systems (human anatomy) is necessary to draw conclusions from a musculoskeletal modelling simulation.
5.4. The role of anthropology in musculoskeletal modelling
Finally, it is important to make clear that in this demonstration we have used a walking trial of one subject and focused on the subject's left side. The kinematics and GRFs produced by this individual's right side and by other walking trials demonstrate a consistent pattern, but were not numerically identical to those reproduced here (electronic supplementary material). Consequently, the muscle forces and joint reactions are also somewhat different. The degree of intraindividual variability in gait parameters is influenced by such characteristics as age and health, but, even for this healthy adult male, peak forces within a muscle can vary. In muscles that exert relatively small forces (less than 500 N), step-to-step values can vary by more than a factor of 2 (electronic supplementary material). Muscles which show greater force generations (e.g. gastrocnemius; greater than 1000 N), step-to-step variability is lower (a factor of approx. 1–1.4; electronic supplementary material). This serves to underscore the importance of collecting and analysing multiple steps when the intent is to report patterns and forces (not the intent here).
While scaling addresses some aspects of variability and improves model predictions [35,69], scaling itself does not solve all the problems. Interindividual variability in such important variables as muscle parameters has been appreciated for over a decade [89], but the consequences of such variability for musculoskeletal modelling remains unclear. For instance, muscle strength as implemented in the solution that we use for this report depends on three factors: muscle volume, optimal fibre length and pennation angle. The MoCap model, which we used in our demonstration, uses a set of values derived from several sources but is critically informed by Klein Horsman et al. [90] and Carbone et al. [34]. These data are from measurements of cadavers who were 77 and 85 years old at the time of death, respectively. Cadaveric data are widely used in musculoskeletal models (e.g. [53,99,117]), but muscle parameters from cadavers may not represent those from individuals who are not at imminent risk of death. Table 1 includes values from Klein Horsman et al. [34] and Carbone et al. [34] and from a recent analysis of 10 healthy individuals who were evaluated with diffusion tensor imaging [91]. Gastrocnemius and soleus, the principal drivers of forward propulsion, vary in all three factors. Importantly, the complexity of gait mechanics that makes musculoskeletal modelling an attractive solution also makes understanding the effect of this variation impossible without simulating the difference. Larger volume, smaller pennation angle and shorter optimal fibre lengths all act to increase muscle strength, thereby increasing the allocation of force to it by the algorithm. But, the algorithm uses a muscle's strength relative to all other muscles' strengths, making it impossible to predict the impact of muscle parameter variation on muscle forces or joint reactions. Variation in muscle parameters is but one of the many ways that people vary in their musculoskeletal anatomy. Leveraging the expertise of biological anthropologists to assess the impact of human biological diversity on musculoskeletal modelling can only be beneficial for the entire community of researchers.
6. Conclusion
Although motivated by different specific goals, biological anthropologists and engineers (principally biomedical) both pursue questions that revolve around understanding the internal forces of the primate musculoskeletal system. Engineering questions tend to be species-specific, focusing on living humans. Biological anthropologists are also interested in living humans, but frequently ask questions that include greater taxonomic diversity. Research in both fields can be specific to the details of the individual (e.g. patient-specific orthopaedics, particular hominin fossils) or generic representatives of groups (e.g. sex-specific running shoes, comparisons of primate taxa).
Musculoskeletal modelling is powerful because of the flexibility it provides at every stage of research. Building a new model requires intensive effort, but may be necessary based on the question. Previously validated models from repositories lower barriers to modelling experiments, but researchers should be aware of choices made by model creators. Like all modelling exercises, serious and intentional thought must be dedicated to the multitude of modelling choices. This is not to say that all models need to be complex. If a simpler model adequately describes the phenomenon of interest and answers the question at hand, then a simpler model may be preferable.
Critical to support rigorous simulations, musculoskeletal models include changeable parameters that describe all aspects of the musculoskeletal system. This is a purposeful feature included by engineers that drove musculoskeletal modelling and an attractive feature for biological anthropologists. Model users can change the size and shape of the skeletal elements, dimensions and inertial properties of body segments, DOFs and range of motion of joints, action of passive soft tissue structures, the geometry and force capacity of muscles, as well as numerous other parameters. Generic models from repositories (such as the one used here) can be scaled and reshaped to capture many aspects of human variation (e.g. limb dimensions, mass), and hence are well suited to many questions. More effort can be dedicated to refine aspects of a model that are of particular interest (e.g. change the morphology of the bony pelvis). Modelling flexibility naturally extends to the information that can be derived from simulations. Forces generated by individual muscles (or portions of muscles) as well as internal to joints are standard output, as is the ability to express these quantities in different coordinate systems. Until internal forces of the body can be measured non-invasively, musculoskeletal modelling will remain the pre-eminent approach to these questions.
Future collaborations between engineers and biological anthropologists to create better models of human movement have the potential to expand the understanding of mobility in both fields. Subject-specific models are internally consistent but not generally representative of people, so extrapolating should only be done with this limitation in mind. Developing a model from scratch for a fundamentally different organism (e.g. a chimpanzee) is a significant investment.
Data accessibility
Additional information is available as the electronic supplementary material. Raw data and simulation output are available at the assigned figshare link.
Authors' contributions
A.D.S. and P.A.K. conceived the project. S.G.L. and P.A.K. designed and carried out the data collection. A.D.S. and P.A.K. carried out the simulations and data analysis. A.D.S. and P.A.K. wrote the draft and all authors contributed to the final version of the manuscript.
Competing interests
We declare we have no competing interests.
Funding
P.A.K. and S.G.L. thank S.K. Benirschke, MD, Jerome Debs Endowment Chair in Orthopaedic Traumatology at the University of Washington, for the ongoing support of this research and for our many conversations about foot form and function.
Acknowledgements
The authors would like to thank the participants of this study for their time and patience during the protocol. The authors would also like to thank the reviewers and editors for comments that improved this manuscript.
Footnotes
One contribution of 12 to a theme issue ‘Biological anthroengineering’.
Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5486386.
Review articles
A review of musculoskeletal modelling of human locomotion
,
and
Published:13 August 2021https://doi.org/10.1098/rsfs.2020.0060
Abstract
Locomotion through the environment is important because movement provides access to key resources, including food, shelter and mates. Central to many locomotion-focused questions is the need to understand internal forces, particularly muscle forces and joint reactions. Musculoskeletal modelling, which typically harnesses the power of inverse dynamics, unites experimental data that are collected on living subjects with virtual models of their morphology. The inputs required for producing good musculoskeletal models include body geometry, muscle parameters, motion variables and ground reaction forces. This methodological approach is critically informed by both biological anthropology, with its focus on variation in human form and function, and mechanical engineering, with a focus on the application of Newtonian mechanics to current problems. Here, we demonstrate the application of a musculoskeletal modelling approach to human walking using the data of a single male subject. Furthermore, we discuss the decisions required to build the model, including how to customize the musculoskeletal model, and suggest cautions that both biological anthropologists and engineers who are interested in this topic should consider.
1. Introduction
A principal research focus within biological anthropology is interpreting skeletal variation within the context of behavioural diversity, including variation in diet, disease processes and activity patterns among living and extinct primates, including humans and hominins [1–4]. Locomotor behaviour and musculoskeletal morphology have been a central focus of this area of inquiry because locomotion provides primates with access to key evolutionary resources including food, shelter and mates. Furthermore, bipedalism is considered a defining character of the hominin lineage, cementing its importance within the field [5–7]. The challenge of elucidating form–behaviour relationships are the complexities of the intervening and underlying functions (e.g. joint range of motion) and biomechanics (e.g. joint reaction forces).
Early anthropological research connecting form and function was both typological and qualitative in nature; however, anthropologists have increasingly adopted mechanical engineering approaches to answer questions about primate locomotion generally, and hominin locomotion in particular (see [4] for a historical review). Within anthropology, two major subfields of mechanics have been applied to locomotor systems: Newtonian mechanics and solid continuum mechanics. Newtonian mechanics is concerned with the description of the motion of solid bodies under the influence of a system of forces, while solid continuum mechanics quantifies the behaviour (e.g. deformation) of solid materials when subjected to forces [8,9]. Lovejoy and colleagues [10–13] were early adopters within anthropology of both Newtonian mechanics and solid continuum mechanics [4]. Newtonian mechanics is routinely used by anthropologists to investigate a variety of questions (e.g. [14–19]). Beam theory, an application of solid continuum mechanics, has been used to estimate the structural capacity of long bones and is regularly exploited by anthropologists (see [4] for a review) to reconstruct hominin locomotor adaptations (e.g. [13–17]) as well as modern human mobility patterns (e.g. [18–24]).
The application of Newtonian and continuum mechanics is, of course, a core competency of mechanical and civil structural engineering, and the methods that anthropologists use derive directly from these engineering fields. The techniques and theories that underlie their application were developed, however, to design and analyse human-made products using inert materials. While such structures and machines tend to be faithfully produced based on design specifications, humans are not clones. Rather, each individual is completed from a standard template that has, as an inherent feature, the ability to respond to the environment. This means that whereas human-made products generally have low tolerance for variance, among humans, variation is the rule. Furthermore, the ‘standard’ template itself varies among ancestral lineages, populations and species because of the evolutionary history of each (e.g. [25–27]). The biomedical engineering world has responded to this with the use of patient-specific models because, in patients, a direct relationship exists between a particular morphology (i.e. that of the patient) and the question to be answered (e.g. which implant configuration limits strain on the bone?) with the musculoskeletal model (e.g. [28–35]). Not all questions, however, concern specific individuals. In anthropology, the focus may be comparing species (or populations) to understand responses to selective pressures. In engineering, running shoes, for example, are not designed for particular individuals, who vary in foot function [36], but rather for people generally. And people are not just scaled versions of each other [37]. Thus, an appreciation of human variation, a core competency of biological anthropology, is a critical component of any modelling exercise.
At the heart of many locomotor-focused questions, both engineering and anthropological, is a need to evaluate internal (e.g. muscle, joint contact) forces. Finite-element analysis of skeletal elements [38,39], estimating locomotor costs [40,41] and evaluating task performance [42,43] all require muscle forces to be known or estimated. Measuring muscle and joint contact forces in living animals is, however, exceptionally challenging because it requires surgical intervention to implant the relevant sensors into the body [44,45]. As an additional obstacle, these techniques are normally limited to a single joint or muscle/tendon within a single study (see [44] for examples). Consequently, internal forces must be estimated. Poorly estimated loads are likely to produce misleading results, interpretations and conclusions. Consequently, estimating internal forces (e.g. muscle forces) in humans is of critical interest to both anthropologists and engineers.
1.1. Musculoskeletal modelling: a path forwards
Originally driven by a clinical interest in gait pathologies, musculoskeletal modelling has emerged as the pre-eminent technique to elucidate the internal details of human (and other animal) movement [46]. These models build upon dynamic analyses of linked rigid-segment models ([47] and described below). Rigid-segment models represent the body as a series of solid body segments which are linked together at joints. Musculoskeletal modelling extends rigid-segment models by including detailed models of individual muscles, as well as models of neurological control. In extending the dynamic approach, musculoskeletal modelling solves the muscle redundancy problem based on muscle parameters under a specified minimization criterion [48–51]. The muscle redundancy problem refers to the fact that there are more muscles than necessary to actuate movements at joints [48,49]. In solving for muscle forces, this leads to an indeterminate system of equations in which there are more unknowns (muscle forces) than constraints (equations) [51], thus requiring optimization solutions. Results of musculoskeletal modelling can reveal individual sequences of muscle excitation and activation, as well as estimates of muscle forces and bone-to-bone contact forces.
In a series of papers, Pedersen et al. (see [45]) developed an early version of this type of model and used it to estimate lower limb muscle forces and hip joint contact forces. The hip joint contact forces were validated using an instrumented hip joint replacement in a single elderly subject [45]. Consequently, these data provided an early demonstration of the validity of musculoskeletal modelling as an approach to determining internal forces. Since then musculoskeletal modelling has been used for numerous projects exploring human motion and locomotion [31,43,52–56] as well as that of non-human primates [57–59] and early hominins [60]. Wang et al. [61] applied this technique to the hominin fossil record, but their approach was necessarily limited by the lack of software available to carry out the complex calculations. Fortunately, several software packages (e.g. AnyBody Modeling System™, OpenSim, FreeBody) are now available to carry out musculoskeletal modelling [45,62,63].
Accurate musculoskeletal models demand both understanding of the principles of Newtonian mechanics and an astute appreciation of anatomy and anatomical variation. Thus, successful musculoskeletal modelling is clearly situated at the interface between engineering (i.e. mechanics) and biological anthropology (i.e. human variation). Here, we provide an overview of this effort from an anthroengineering perspective, blending knowledge from both domains with a special emphasis on a task—understanding human walking—that is of equal interest to both groups.
2. Background
2.1. The model
In musculoskeletal modelling used to solve for muscle forces, the body of interest (e.g. the whole human body) is modelled as a series of rigid segments connected by joints and actuated by muscles [34,46] (figure 1). This body moves in a virtual space whose dimensions are defined by a global Cartesian coordinate system. Models can include as few or as many body segments as necessary to answer the question of interest. For walking, models minimally include segments for the pelvis (and trunk/arms/head, often aggregated) as well as left and right feet, shanks and thighs. More complete models often also include segments that represent the patella and talus, as well as a head–arm–trunk (HAT) segment separate from the pelvic segment. Additionally, the upper part of the body may also be modelled to include segments for the lumbar region, thorax, arm, forearm, hand, neck and head segments. Body segments are defined in terms of their physical dimensions, mass, location of centre of mass and moments of inertia. Each body segment also includes its own Cartesian coordinate system, usually aligned with standard anatomical directions of the segment (e.g. proximodistal, anteroposterior, mediolateral). Other coordinate systems can be defined relative to that of the segment or to the global system and allow for ease in defining model geometry. The segments further rely on the assumption that they do not deform during locomotion (i.e. are rigid). In the software modelling/simulation environment, body segments are almost universally represented by virtual surface models of the relevant skeletal element(s).
Figure 1. MoCap model mid-simulation of walking. The image shows the model (rigid body segments represented by skeletal elements and visual representation of muscle models), four force plates (grey), the ground reaction force vector (blue line through model), force plate readings (blue lines below the force plate), experimental marker locations (blue spheres), virtual markers (red spheres with coordinate system arrows) and global coordinate system (yellow arrows).
Body segments are connected to other segments at joints, which represent articulations (e.g. hip, knee, ankle) within the body. Joints are frictionless and non-deformable and, consequently, forces and moments transfer across these joints in an equal and opposite manner. A joint definition includes its location relative to the body segments it connects, as well as the degrees of freedom (DOFs) permissible at the joint (one, two or three axes of rotations and, less commonly, translations). The direction of rotation axes are usually expressed in terms of a local (body segment) coordinate system.
With the body segments and joints defined, the skeletal portion of the musculoskeletal model is sufficient to conduct dynamic analyses. Musculoskeletal models extend the analysis by including models of individual muscles or even portions of muscles (figure 1). Muscle models range in their complexity and, as with other parts of the model, the necessary complexity depends on the motion being simulated. (For a description of how muscle complexity can be modelled, see §5.3.) A muscle model minimally requires geometry (i.e. path definition) and a parameter representing its maximum possible force. More complex muscle models also include parameters for pennation angle, optimal fibre length, physiological cross-sectional area (PCSA) and activation and deactivation times. Muscle geometry must include the coordinates of its origin (location relative to the proximal segment) and insertion (location relative to the distal segment), but muscles can be geometrically more complex than a straight line path between origin and insertion. The most basic of these complexities derives from a muscle with a broad origin (e.g. gluteus medius) or insertion (e.g. adductor magnus). In these muscles, a single origin-to-insertion path does not encompass their possible actions. One way to solve this is to split the muscle into muscle parts, referred to here as muscle elements. Each muscle element represents a region of the muscle (e.g. the superior portion of adductor magnus) and can have its own set of muscle parameters.
Most modelling systems allow for additional coordinate positions that are used to define the curved path that the physical muscles take between bony attachments (e.g. the path sartorius takes across the anterior thigh as it traverses from anterior superior iliac spine (ASIS) to proximal tibia). Via points, the simplest type of path point, are one or more coordinate positions that are intermediate between the origins and insertion points. The muscle assumes a series of straight line paths between the origin, via point(s) and insertion. The location of the via point is fixed relative to the segment's centre of mass and the connection to other segments. With via points, the muscle is required to pass through each specific point at all times. These points allow muscles to pass over other body surfaces (e.g. iliacus wrapping over the pelvic brim). In addition, muscle geometry paths can be altered by the use of ‘wrapping surfaces’ that are designed to account for relative motion between musculoskeletal elements and other soft tissue (e.g. rectus femoris sliding across the soft tissue structures anterior to the hip). These wrapping surfaces are geometric objects (e.g. cylinders, ellipsoids) that are fixed relative to a particular body segment. Muscle paths are not constrained to pass through any particular coordinate point (as with via points), but rather can move freely over the wrapping surface to find the shortest possible path across the surface between origin, insertion or via points that are beyond the wrapping surface. As the proximal and distal segments change orientation relative to each other, the point of contact between the muscle and the adjacent musculoskeletal elements can move (e.g. gastrocnemius relative to the posterodistal femur as the knee flexes). Forces developed from the action of the muscle act on the body segment used to define the geometry (origin, insertion, via points and/or points on wrapping surfaces). Muscles may also be modelled to be more complex in terms of the number of parameters that define force generation. Minimally, a muscle model must include a maximum force value which defines the force when the muscle is maximally activated, but other approximations of muscle force generation are available, e.g. the Hill-type muscle model [64].
The final required model element is a set of virtual motion markers, which is needed to match model motion to motion measured experimentally (figure 1). For motion tracking (described below), the human subject is fitted with a set of markers that are tracked by a camera system. The musculoskeletal model must include one virtual marker for each experimentally measured marker. Each virtual marker is connected to a body segment and defined in terms of that segment's coordinate system. For instance, if researchers affix a motion-tracking marker over the ASIS of a human participant, the model will require a virtual marker connected to the pelvis body segment with coordinates that position it correctly relative to the model ASIS.
2.2. The input data: motion data and ground reaction forces
An inverse dynamics approach to musculoskeletal modelling (described below) requires that the motions (translation and rotation) of the rigid segments (i.e. the skeletal elements) and externally applied forces be known. This process is well established and has been described in detail elsewhere (e.g. [65–67]). To obtain motions, retroreflective (infrared) markers attached to the skin or clothing are frequently used. Each body segment requires a minimum of three markers that are trackable throughout the entire sequence of movements, so this traceability should be checked in all portions of the data collection volume. The three-dimensional ground reaction force (GRF) for a complete stance is also necessary. GRFs are typically acquired through force plates embedded either in the floor or in a treadmill.
2.3. The simulation
2.3.1. Scaling and marker optimization
Model scaling and virtual marker position optimization are the processes of determining the model segment dimensions and virtual marker positions so that the virtual model matches the human subject and experimental protocol (motion marker placement on the human subject) as well as possible [68]. The degree of matching between the model and experimental data is evaluated using a global error metric between experimental and virtual motion marker positions. During the inverse kinematic analysis (described below), musculoskeletal modelling software finds the position of the model, including overall location in space and individual joint positions, that minimizes the sum of squared deviations between the virtual marker set and experimental motion marker data. The goal of marker optimization and body scaling is to minimize this error term.
The human models that are included with many modelling software packages are generic human models. For analysis, the model must match the subject for whom the input data are known in terms of mass, limb parameters (e.g. distances between joint centres and segment moments of inertia) and virtual marker positions [69,70]. Marker position optimization and scaling are intertwined processes because marker positions can, in part, determine segment dimensions. Initial scaling of the generic virtual model achieves a loose match between the overall bodily dimensions of the subject and those of the virtual model. If a virtual model from a repository is used (rather than creating a virtual model from scratch), then the model will have been developed using the parameters of the individual whose data were available (i.e. a generic model). This generic model will almost certainly not represent the dimensions of the subject on whom motion-tracking data and GRFs were collected, so the generic model needs to be scaled [69]. The simplest scaling approach is uniform scaling of lengths and volumes. Total body mass is the most basic dimension needed, and the ratio of the body mass of the subject to that of the generic model is typically used to adjust segment masses and other parameters that are volume-dependent. The ratio of subject to generic model stature, a length scale, is used as an initial estimation of segment dimensions and other dimensions that are length-dependent. Other anthropometric data (e.g. distances between landmarks) can be used instead of, or in addition to, stature. Other scaling methods are also available, including fully manual scaling, using body fat to adjust mass distribution among the segments, and hybrid scaling (e.g. manually scaling some parameters and algorithmically scaling others). Furthermore, depending on the sophistication of the software used to complete the mechanical analysis, these initial scalings may be adjusted during marker optimization.
Different software packages conduct marker optimization and secondary scaling of segments using different algorithms and require varying levels of interaction/input from the user. At one end of the user-input spectrum, some systems require that the user optimizes the locations of the markers on the model ‘by hand’, moving them relative to their parent segment and then interrogating error values associated with each marker during motion tracking. This process is tedious and may require significant effort before optimal marker position is achieved. At the other end of the spectrum, some software packages change both segment dimensions and virtual motion-tracking markers automatically in a single cohesive optimization process that requires little user interaction. In these frameworks, the user can permit or restrict changes in virtual marker position in none, some or all of the coordinate directions. For instance, motion-tracking markers associated with easily palpable bony landmarks (e.g. medial malleolus) are also easy to locate on the virtual model (referencing the virtual skeleton) and may be considered ‘fixed’ in their location relative to their body segment (the shank). Thus, during the optimization process, such a fixed marker would not change its position relative to the segment, but would influence the determination of shank length. Other markers (e.g. marker plates attached mid-segment) may be allowed to optimize their position in all three coordinate directions, and then contribute little to segment length determination. Finally, some markers may be well defined in two coordinate directions, but less so in the third. For instance, the superoinferior and mediolateral positions of an ASIS marker may be well defined and constrained from moving relative to the pelvis during the optimization process. The anteroposterior position of the virtual ASIS marker relative to the virtual pelvis may be more difficult to estimate as it reflects tissue thickness over this bony landmark (which is hard to measure on living subjects in the absence of internal imaging modalities). Thus, the user might fix the ASIS marker in two coordinate directions, but allow its position to optimize in the third. The fixed mediolateral and superoinferior position on the two ASISs would then influence the scaling of pelvic width.
2.3.2. Inverse kinematics/motion tracking
To carry out inverse dynamics and, in turn, estimate internal forces and moments, segment linear and rotational accelerations must be obtained. Most musculoskeletal modelling software (e.g. Anybody Modeling System™, OpenSim) uses a kinematic analysis procedure of over-determined biomechanical systems, which is formulated as a weighted least-squares optimization problem that matches experimental and model marker positions (e.g. [47,71]).
2.3.3. Inverse dynamics
While statics is the branch of mechanics that deals with stationary objects, dynamics is the branch in which motions are the focus. In statics, an equilibrium that greatly facilitates examination of the system exists where the environment (e.g. the substrate) and the structure (e.g. the human body) are balanced. In standing, the body pushes through the feet and against the ground exactly as much as the ground pushes against the body through the feet. In dynamics, the balance is time-dependent and the solution to the system is through the equations of motion, which relate forces to motions. When the forces that produce the motions are fully known, the problem is characterized as forward dynamics, while when motions are known it is described as inverse dynamics.
In human bodies, the measurement of motions (via marker tracking or other means) is much less invasive and more feasible than measuring the muscle forces that produce the motions, so currently most analyses that depend on the musculoskeletal forces generated in human gait are typically treated as inverse dynamic problems (e.g. [72–74]). Further complicating the solution in multi-body problems (i.e. those with multiple, linked rigid bodies) is that the system can be interrogated from an external or internal perspective. The external view includes the action of the GRF to produce the movement of the body along a path, while the internal one includes both the joint reactions and the action of muscle forces in creating changes in joint angles. In mechanics, the relationship between the number of unknown variables and the number of known variables determines the difficulty of the solution. Fortunately, the internal and external perspectives can be exploited to develop solutions. For instance, in the single-stance foot, the GRF is the external force and produces the moment acting on the foot segment and reacted to by the ankle joint reaction force and moment (which are an internal force and moment). Once the ankle joint reaction force and moment are known, they can be combined with the motion of the shank to solve for the knee joint reaction and moment. The inverse dynamic solution allows for this time-dependent chain solution to determine the joint reactions.
2.3.4. Muscle redundancy solution
A system where the number of variables with unknown value equals the number of independent equations that fully define the system is determinant and a unique solution exists. For systems where there are more variables with unknown values than equations—called indeterminate—assumptions are made to provide more equations [75,76]. Many biological problems, including developing muscle forces, are indeterminate because there are more muscles capable of acting across a joint than necessary from a purely mechanical perspective, i.e. some muscles are (mechanically) redundant. If the desired result from an analysis that uses inverse dynamics is the development of muscle forces, then the joint moments determined from the inverse dynamic solution need to be apportioned to individual muscle forces. This is not a trivial undertaking, because the degree of redundancy is often quite large, so no general consensus on the most accurate criterion on which to base an apportionment algorithm currently exists [77].
Musculoskeletal modelling solutions typically include an algorithm that uses muscle parameters to apportion the joint moments to particular muscles (as muscle forces). Early solutions used PCSA to solve the muscle redundancy problem [75]: the apportionment was made based on the ratio of the PCSA of a particular muscle relative to the total PCSA of all muscles capable of acting (e.g. [61]). While effective from an algorithmic perspective, PCSA-based apportionment does not reflect the natural solution, because cross-sectional area approaches assume equal activation of all possible muscles.
Another algorithm for solving the muscle redundancy problem attempts to minimize the sum of muscle activations raised to a power [76]. The activation represents a proportion of total maximum strength the muscle can generate given the current state of the model (i.e. activation is generally bound between 0 and 1). If the muscle activation exponent is 1, i.e. a linear combination, force is allocated to the muscle with the greatest moment arm until it reaches its maximum allowed value, and then other muscles are recruited. Exponents greater than 1 will distribute forces more evenly across muscles available to produce the moment at a particular joint. The greater the power, the more even will be the distribution (PCSA-based activations are essentially an exponent of infinity). Other criteria for apportionment of the joint moments have also been advanced, such as electromyography (e.g. [78,79]), but activation-based methods remain the standard [76]. Until in vivo muscle forces can be reliably measured, all solutions to the muscle redundancy problem should be treated as hypotheses about how joints and muscles produce motion.
2.4. The output
A valuable aspect of musculoskeletal modelling (or of any model) is that variables of interest can be interrogated anywhere within the model. Standard output from a musculoskeletal model includes linear and angular displacements, velocities and accelerations for all body segments. Similarly, joint positions, velocities and accelerations are also tracked and can be exported for analysis. Muscle length, activation, moment arm and force are all possible outputs, as are joint reaction forces. Additionally, users can track the position of muscle origins/via points/insertions throughout the modelled motion. Conveniently, joint reaction forces can be output in the global coordinate system or the coordinate system of the relevant body segments.
Human locomotion requires the repetition of cycles of limb and body motions that can be reduced to a stride. For walking, a step is defined from an event (e.g. initial contact) on one foot to that same event on the contralateral foot and a stride is two contiguous steps. A stride can vary in duration, even at a constant velocity, so creating an average stride or evaluating the variability of stride parameters requires normalization of the strides to some consistent duration. The typical choice is per cent of stride cycle, where the stride duration is described in terms of equal portions of the entire stride (e.g. [80,81]). Evaluating the variability of parameters of interest (e.g. muscle force, joint reactions) can be an important step in validating the calibration of a solution. While some differences among trials of the same individual are expected, an evaluation of the variability in light of the question that the model is designed to answer should occur.
2.5. Human variability
A crucial advantage of musculoskeletal modelling is the opportunity to include aspects of morphological variation, including the production of subject-specific models. Model scaling and motion tracking incorporate some aspects of human variability, including differences in body mass, limb lengths/dimensions and kinematics, but other aspects are less obviously accessible for modification in musculoskeletal modelling software. These less accessible aspects, such as variation in bone or muscle morphology, are also important to consider. Quantifying human morphological (musculoskeletal) variation is at the heart of biological anthropology. For example, human (and hominin) pelvic morphology has a profound influence on bipedal capabilities [14,80], and has been the focus of intense scrutiny [41,81,82]. Humans are known to be sexually dimorphic in pelvic size and shape [83–85]. Furthermore, human pelvic morphology is also known to vary by population [86–88]. Sources of variation can be incorporated into musculoskeletal modelling to understand pathology [31] and differences in joint reaction forces [35]. For example, marker-informed parameter optimization and secondary scaling (described above) can change the width of the pelvis if experimental and virtual markers are included on bony landmarks of the pelvis, such as the ASISs. If the subject is wider than the generic model, the parameterization can widen the pelvis so that the virtual and experimental markers are more closely aligned. Scaling can also be accomplished manually either by applying factors to the three segment coordinate directions or by redefining the locations of specific anatomical landmarks associated with the segment. The former method moves all points proportionately while the latter allows for localized warping.
In addition to skeletal variation, aspects of muscle morphology that vary among humans may also be critical for accurate musculoskeletal modelling. For instance, variability in muscle parameters has been of concern for some time [89]. While repository models represent a particular person's morphology, people differ in important musculoskeletal modelling parameters such as muscle volume, optimal fibre length and pennation angle (table 1).
Table 1. Comparison of muscle parameters between the standard MoCap model in AMMR™ v. 2.3 (AnyBody Technology, Aalborg, Denmark; from Klein Horsman et al. [90] and Carbone et al. [34]) and Charles et al. [91]. Muscle names are per the MoCap AMMR v. 2.3. Where multiple muscle regions are present in the MoCap model, the values for the muscle volumes are summed while the optimal fibre lengths and pennation angles are averaged.Collapsemuscleoptimal fibre length (cm)pennation angle (°)muscle volume (ml)MoCapCharles—averageMoCapCharles—averageMoCapCharles—average
adductor brevis | 10.38 | 7.6 | 0.0 | 11 | 108.90 | 93 |
adductor longus | 10.57 | 11.0 | 0.0 | 12 | 159.56 | 159 |
adductor magnus | 9.97 | 23.1 | 0.0 | 12 | 559.25 | 567 |
biceps femoris caput breve | 9.14 | 10.9 | 0.0 | 9 | 107.95 | 92 |
biceps femoris caput longum | 8.54 | 20.4 | 29.9 | 11 | 232.01 | 194 |
extensor digitorum longus | 5.99 | 13.8 | 8.3 | 7 | 32.29 | 76 |
extensor hallucis longus | 5.98 | 10.6 | 14.4 | 7 | 36.27 | 21 |
flexor digitorum longus | 3.84 | 28.4 | 25.28 | |||
flexor hallucis longus | 5.20 | 30.1 | 161.50 | |||
gastrocnemius lateralis | 5.69 | 12.2 | 25.4 | 9 | 136.36 | 128 |
gastrocnemius medialis | 6.01 | 9.7 | 10.8 | 10 | 263.26 | 230 |
gemellus | 3.43 | 0.0 | 28.40 | |||
gluteus maximus | 13.57 | 0.0 | 936.55 | |||
gluteus medius | 6.76 | 5.3 | 674.71 | |||
gluteus minimus | 3.28 | 0.0 | 82.59 | |||
gracilis | 18.11 | 17.3 | 0.0 | 6 | 87.97 | 91 |
iliacus | 8.16 | 8.8 | 203.13 | |||
obturator externus inferior | 6.88 | 0.0 | 37.88 | |||
obturator externus superior | 2.77 | 0.0 | 68.18 | |||
obturator internus | 2.05 | 0.0 | 52.08 | |||
pectineus | 9.50 | 0.0 | 64.58 | |||
peroneus brevis | 4.54 | 23.1 | 86.09 | |||
peroneus longus | 5.08 | 15.8 | 121.52 | |||
peroneus tertius | 4.27 | 19.1 | 26.52 | |||
piriformis | 3.88 | 0.0 | 31.25 | |||
plantaris | 4.77 | 0.0 | 11.36 | |||
popliteus | 2.40 | 7.4 | 0.0 | 8 | 25.57 | 15 |
psoas major | 9.92 | 13.4 | 193.18 | |||
psoas minor | 5.91 | 0.0 | 6.63 | |||
quadratus femoris | 3.37 | 0.0 | 49.24 | |||
rectus femoris | 7.84 | 14.2 | 21.9 | 8 | 226.33 | 249 |
sartorius | 34.71 | 40.8 | 0.0 | 205.49 | 143 | |
semimembranosus | 8.09 | 15.8 | 25.0 | 12 | 138.26 | 244 |
semitendinosus | 14.16 | 18.3 | 0.0 | 8 | 208.33 | 186 |
soleus | 4.40 | 14.6 | 61.6 | 12 | 792.91 | 461 |
tensor fasciae latae | 9.49 | 0.0 | 83.33 | |||
tibialis anterior | 6.83 | 13.7 | 9.6 | 7 | 181.49 | 129 |
tibialis posterior | 3.78 | 34.2 | 163.36 | |||
vastus intermedius | 7.68 | 18.1 | 11.8 | 11 | 292.61 | 521 |
vastus lateralis | 6.69 | 19.6 | 0.0 | 15 | 583.33 | 606 |
vastus medialis | 7.16 | 15.9 | 0.0 | 14 | 454.27 | 415 |
Several possibilities for future work that explores the impact of morphological variation on muscle and joint forces are worth noting. For instance, while some anthropological work has elucidated the interactions of morphology and walking conditions (e.g. velocity, changes of direction, burdens and gradients) on energy expenditure (e.g. [92–94]), much less is known about how morphology and walking condition factors interact to affect joint and muscle forces (although see [31] for an example). People vary morphologically and walking under varied conditions is an activity of daily living, so this area is an important one for many disciplines to pursue including biological anthropology, biomedical engineering and product design. Perhaps more importantly, collecting the requisite data to inform a musculoskeletal model is time consuming, invasive and dependent on (at this time) fragile laboratory equipment. Consequently, acquiring samples that represent the breadth of human morphological variation is unlikely to occur in the near future. Once validated, though, a musculoskeletal model can be modified to assess the impact of these known morphological differences on muscle and joint forces. Extending that line of research further, a longer-term goal would be use musculoskeletal modelling to understand the influence of the morphological differences between modern humans and extinct hominins on muscle and joint forces, an area that has sought biomechanical solutions without effective software support for decades [12,14,61,95,96].
2.6. Current study
Here, we demonstrate the process and output of one human subject during walking using AnyBody Technology's musculoskeletal modelling software (Anybody Modeling System™; AnyBody Technology, Aalborg, Denmark). We present several of the possible outputs that could be of interest to biological anthropologists. This includes joint kinematics, joint reaction forces and a selection of lower limb muscle force magnitudes. We make some recommendations for future researchers to help facilitate their use of software and output. Additionally, by presenting a set of possible outputs, we demonstrate the motivation for undertaking musculoskeletal modelling and that results may have a direct bearing on anthropological questions, or as intermediate steps to other analyses (e.g. finite-element analysis). Additionally, while engineering has much to offer biological anthropologists, the latter have an appreciation for human variation which may enhance future generations of musculoskeletal models.
3. Material and methods
3.1. Participant
One healthy participant was used for this demonstration, but he is part of a larger study aimed at examining human walking. The participant was 36 years old (male; body mass, 74.3 kg; stature, 172 cm) and reported being free from lower-limb injuries. The University of Washington's Institutional Review Board approved all aspects of this study (IRB no. STUDY00001125).
3.2. Experimental protocol
Kinetic and kinematic data were measured using a four-force plate (Kistler, Switzerland), 10-camera motion capture system (Qualisys, Sweden) in the Amplifying Movement & Performance (AMP) laboratory of the University of Washington. Thirty infrared-retroreflective 5 mm hemispherical markers were affixed to anatomical landmarks and used to track motion through the laboratory (figure 1; electronic supplementary material). The participant walked unshod the length (10 m) of the gait laboratory five times at his self-selected preferred walking speed (1.15 m s−1). The participant walked in a straight path with his eyes directed at a point on the far wall of the laboratory and contacted the surface of each force plate with a single foot. Trials with multiple feet on a single force plate or ones where the stance foot crossed plates were immediately redone. The participant was allowed several attempts prior to data collection to become familiar with the protocol. All trials exhibited similar kinematics and kinetics (data not shown), so we selected one trial to use in this demonstration.
Marker data and GRFs were collected at 120 Hz and 1200 Hz, respectively, and then were filtered at 5 Hz and 10 Hz, respectively, with a fourth-order, low-pass zero-lag Butterworth filter. Calibration of the system yielded a limitation in its fidelity for marker data of 1 mm and force data of ±2.5 N for the direction of travel (X), ±5 N side (Y) and ±25 N vertical (Z). All data were exported from the Qualisys Track Manager software in C3D file format, which can be read directly by the modelling software (provided in electronic supplementary material).
3.3. Musculoskeletal model
We used the ‘MoCap’ model from the AnyBody Managed Model Repository™ (AMMR v. 2.3; AnyBody Technology, Aalborg, Denmark; http://www.anybodytech.com/software/model-repository-ammr), which is a multi-trial, full-body, motion-capture-driven human gait model, and the commercially available AnyBody Modelling System™ (v. 7.3; AnyBody Technology, Aalborg, Denmark) to calculate forces in the lower limb [34,97]. The MoCap model that we used for this demonstration has been used to assess human gait in numerous applications (e.g. [98,99]). The MoCap model is assembled from a trunk and left and right lower limb components. The trunk includes a lumbar region, trunk, neck and head. The lower limbs consist of the following segments: thigh, patella, shank, talus and foot. Each lower limb has six total DOFs, including all three rotations at the hip and one each at the knee (flexion/extension), ankle (plantar flexion/dorsiflexion) and subtalar (inversion/eversion) joint. The pelvis (relative to the ground) has six DOFs, three translational and three rotational. Forty-one anatomically defined, bilateral muscles (table 1) are represented by 169 muscle elements in each model lower limb (e.g. gluteus medius is composed of 12 separate muscle element actuators) [34]. The muscle elements in this model generate force as a function of their activation and a strength parameter which is appropriate for human walking [97,100]. We defined a virtual motion-tracking marker for each of the experimentally measured motion markers (electronic supplementary material). The polynomial algorithm that was employed to solve the muscle redundancy problem apportioned force to muscles using an objective function with power = 3. (For more information regarding the model, see §5.3.)
3.4. Musculoskeletal simulation
The first step of the simulation was to scale the model to match the dimensions of the subject while simultaneously optimizing the marker coordinate and joint rotation axes' locations [68]. After initial scaling, the values of these parameters (i.e. segment dimensions, marker locations, rotation axes' locations) were determined using the marker motion data from a trial. Values for all parameters were optimized to minimize the sum of square distances between model marker positions and experimental marker positions. Following convergence of parameter optimization, we conducted an inverse kinematics analysis. Finally, we conducted the inverse dynamic analysis, which includes apportioning muscle forces using the algorithm that minimizes the cost function as a sum of muscle activation values raised to the third power. This value has been shown to allocate muscle activations in a way that reflects experimental data [77,100,101].
3.5. Output variables and data processing
We output several kinematic and kinetic variables, including joint kinematics, filtered GRF components and joint reaction forces (electronic supplementary material). Additionally, we output the force magnitude for a selection of the muscle elements for the left lower limb during an entire stride cycle of walking. Muscle element forces within anatomically defined muscles were summed for presentation (e.g. 12 elements that represent gluteus medius). All time-related curves were resampled to 1% increments of the stride cycle.
4. Results
GRF component curves across the gait cycle are presented in figure 2. As is typical for human walking, the vertical component has two peaks (early stance peak, 719 N; late stance peak, 778 N) on either side of a valley (633 N) which occurs around the middle of stance phase. The average of the two peaks is 103% of body weight (729 N), while the valley represents 87% of body weight.
Figure 2. Filtered ground reaction force components. Grey vertical line represents moment of toe-off.
Pelvic rotation and obliquity and hip, knee and ankle kinematics are provided in figure 3. All kinematics are of the left lower limb throughout a single stride. Pelvic motion is described in terms of the right hip (previous stance limb) relative to the left hip (current stance limb) [102]. Hip, knee and ankle joint reaction forces, which include contributions from muscle forces, are provided in figure 4. In each joint, the vertical joint reaction force is greater in magnitude than the anteroposterior and mediolateral components. All three joints show a large peak during the second half of the stance phase, and this late peak is significantly larger in the ankle. Muscle force profiles for 18 anatomically defined muscles across the stride cycle are provided in figure 5. During the stance phase of walking, gastrocnemius, soleus and gluteus medius generate the largest forces. All output data are provided in the electronic supplementary material.
Figure 3. Angular kinematics for the pelvic segment and hip, knee and ankle joints. Vertical grey line represents toe-off.
Figure 4. Joint reaction forces expressed in the global coordinate system. Vertical grey line represents toe-off.
Figure 5. Muscle forces for 18 anatomically defined muscles.
5. Discussion
5.1. Model output
The subject that we used to demonstrate the musculoskeletal modelling approach to modelling walking is an adult male without gait pathology and within species-typical anthropometric characteristics. Consequently, the model-derived motions should accord with those determined experimentally and extensively described in the last 40 years of gait research (e.g. [66,103–105]). The model results presented above do meet this requirement, as discussed below in detail. We focus on the model-derived kinematics, because these kinematics are the foundation for all subsequent analyses in musculoskeletal models. Any analysis should begin by verifying that model-derived variables align with those obtained empirically, and the cause of any discrepancies should be understood or noted as a potential limitation. We belabour the details in the following paragraphs in order to demonstrate the basic approach to model validation.
In our subject at heel strike, the pelvis is rotated approximately 10° in the horizontal plane such that the right hip is posterior to the left hip (e.g. the current stance limb). The pelvis then swings through a neutral position and, while leading up to the next right heel strike, into an orientation with the right hip anterior to left hip (approx. 10°). The pelvis drops in the frontal plane with the right hip inferior to the left (stance) hip soon after heel strike as the right lower limb moves into the swing phase with pelvic obliquity ranging between ± 4°. The pelvis remains in this position until the subsequent right heel strike, after which the pelvis approaches the neutral position. Model-derived pelvic rotation and obliquity are similar to empirically determined values (e.g. [103]).
In the sagittal plane, the left hip begins the stance phase in flexion of approximately 17° and moves through a neutral position and into extension of approximately 17°. Winter [66] describes a less symmetrical flexion–extension cycle at natural cadence (mean values of approx. 20° of flexion and 11° of extension), but the model-derived hip rotations are well within one standard deviation of his mean values. During most of the stance phase, the left thigh is adducted less than 6° in accord with others (e.g. [103]).
While the pelvic and hip joints have three rotational DOFs, the knee and the ankle only rotate in the sagittal plane. The left knee flexes slightly (10°) early in the stance phase, but maintains a near fully extended position through midstance, after which it begins to flex in preparation for the swing phase where flexion exceeds 60°. After the initial heel strike and as the forefoot contacts the substrate, the ankle plantar flexes slightly (less than 5°) and then assumes a dorsiflexed position for most of stance. Rapid plantar flexion occurs in late stance as the body is propelled forwards, reaching a plantar flexion value of greater than 10°. The ankle returns to a dorsiflexed position for swing. These values for both knee and ankle angles are typical (e.g. [66,103,105]).
Consequently, the kinematics of the model match those that we expect empirically. While this check is not sufficient to ensure that the internal (muscle and joint reaction forces) are reasonable, it is a necessary condition and validation for any model. Incorrect kinematics guarantee incorrect or suspect muscle forces.
Joint reaction forces have been empirically measured only with instrumented joint implants. By definition, the surgery creates conditions that limit applicability of the results, but no other method currently exists to measure force. Pedersen et al. [45] measured hip joint contact forces in a 72-year-old man (63.2 kg) with a hip replacement and found that the vertical force was 1750 N, considerably less than the almost 3500 N we found with the musculoskeletal model. Important differences between this patient and our subject exist, including age (72 versus 36 years), body mass (62.3 versus 74.3 kg), health status (58 days post hip replacement surgery versus healthy) and walking speed (0.89 versus 1.15 m s−1). While hip and knee joint reaction forces presented here are larger than in vivo values using instrumented prostheses [77,106,107], our model-derived joint forces are, however, similar to other models of healthy adults [35,107].
Although it is difficult to compare muscle forces with experimental data, it is possible to consider the pattern of activation with reference to the muscle's action. The muscles that generate the highest force in the model are gluteus medius, soleus and gastrocnemius. Gluteus medius plays a critical role in bipedal walking by preventing pelvic drop during single stance [103,105]. The gluteus medius muscle profile from the model is double peaked, and the timing of the two peaks is coincident with the two peaks in the vertical GRF, indicating the muscle's support of the pelvis during single stance. Gastrocnemius and soleus plantar flex the ankle; when the foot is in contact with the ground and the hip is extended, plantar flexing the ankle moves the body anteriorly [103,105]. The gastrocnemius and soleus muscles show the greatest force generation during the propulsive portion (second half) of the stance phase, in keeping with our expectation from their function.
Other muscles also display moderately high forces. Tibialis anterior, which dorsiflexes the ankle, shows muscle force generation immediately after heel strike when the body begins to move over the foot, and then is mostly inactive until the second half of the swing phase. The quadriceps demonstrate a force peak early in stance phase, and again at the stance–swing transition, potentially relating to their role in extending and stabilizing the knee [105]. The hamstring muscles show the greatest peak at the end of swing phase and into early stance phase. Activation of these muscles at the end of swing phase is associated with decelerating the limb prior to heel strike. Early in stance the GRF is directed anterior to the hip and knee, and the hamstring muscles are active to prevent further flexion (hip) and extension (knee) of these joints. The hip flexors are active late in the stance phase, and may be active to counteract hamstring hip extension moment as the hamstrings flex the knee in preparation for swing phase. The adductors show initial peak in stance (adductor magnus) and late stance (adductor longus), but generally smaller than forces in other muscles.
5.2. Input data choices
The model-derived motions and joint reaction and muscle forces align with our expectations from the experimental data in the literature. We achieved these results through careful attention to detail when we collected the data we used as inputs to the musculoskeletal model. We provide below some brief suggestions from our experience. Because motions are derived from the change in position of these markers, marker placement on the segment should be chosen to define fully the motion of the segment and to reduce marker positional noise as much as possible. Positional noise arises from the inherent error of measurement of the marker's position in space (system error) and from movement of the surface (skin or clothing), where the markers are placed relative to the underlying rigid segment (e.g. [108–110]). System error is reduced through calibration; for our system, it is less than 1 mm for any marker. Movement between the external (skin) surface and the bone is more complicated to reduce, but may not be clinically relevant [110].
Marker position is a compromise with soft tissue coverage of underlying bone, landmark palpability (to aid placement) and distance between markers (such that marker positional error is much less than the distance between markers) all being important. The shank and foot have palpable landmarks that are covered by minimal soft tissue (e.g. the malleoli), but the thigh (e.g. greater trochanter) and, especially, the pelvis are more challenging segments on which to affix markers. This problem is exacerbated when the subjects of interest have more soft tissue (muscle or fat) overlying the bony landmark, and considerable soft tissue coverage is common in the hip–pelvis region. One approach to ameliorate some of this difficulty is to use a marker plate which has four (or more) markers mounted on a hard plastic plate that is then strapped to the segment. If a plate can be affixed firmly to the segment, as it can be to the thigh, the marker plate eliminates between-marker error. Unfortunately, marker plates are difficult to affix firmly to the pelvis. Marker placement issues, consequently, are most difficult to overcome to track the motion of the pelvis. Yet, for musculoskeletal models of human locomotion, pelvic motions are critical components of system performance. We achieved the motions shown in figure 3 through extensive protocol development and testing.
Another critical input to the simulation are the GRFs. It is possible to conduct an inverse dynamics solution of the single support phase of walking with a single force plate (e.g. [111]), but extending the analysis to the double support portions of the stride is more difficult with one force plate because it would require estimation of the unknown contribution of the other stance foot. Consequently, any research that requires information about double support should optimally be conducted with multiple force plates. The placement of the force plates relative to each other should be appropriate to the population and gait. Stride lengths for adults or running gaits are longer than for children or walking, and plates should be positioned to prevent subjects ‘targeting’ foot placement or large numbers of repeated trials. Data analysis is less complicated if stance occurs on one plate, but the data from two immediately adjacent plates can be used as long as there is no gap between them [112]. Multiple feet contacting a single plate is also possible to resolve with additional information (such as localized pressure sensors) or predictive methods (e.g. [113,114]). All of these remediation techniques are time-intensive. The laboratory in which the presented data were collected has four floor-embedded force plates that are positioned to allow normal stride lengths for adults when walking. The subject and trial that we used for this demonstration is part of a much larger study [115]. Nonetheless, we assessed each foot placement on all four plates for each trial before moving on to the next trial. Such data checking is a tedious effort at the time of data collection, but well worth the effort to ensure all data are suitable for future analyses. Using four force plates, we had data frames before and after the data that we used in the analysis. This allowed us to discard the initiation and termination frames of the simulation, which are frequently unusable. Consequently, researchers are advised to locate force plates well and test the placement prior to data collection.
Finally, researchers may consider collecting complete anthropometrics of subjects and using these in the initial scaling as the starting points for the parameterization. People vary considerably in terms of their proportions [35,37], so standard models, such as the model we used here, represent a morphology that may not match the dimensions of any particular subject. If your model starts with segment lengths that are substantially different from the segment lengths of the model, the parameterization algorithm may not be able to converge on a suitable solution. Two failures can occur: the parameterization can fail to converge, which produces an obvious error message, or the parameterization converges, but, in order to converge, the simulation requires changes to the morphology or motions that are inconsistent with the experimental data or subject dimensions. This latter failure can be insidious. An absolutely critical step in model validation is to check the motions and shapes of the segments after parameterization. A poorly parameterized model is incorrect and will not produce acceptable motion or joint and muscle forces.
5.3. Modelling choices
In order to produce this demonstration, we made a number of choices with regard to the model which bear discussion. The most important of these was the choice to use a standard, repository model, in this case the MoCap model of AMMR v. 2.3. By using the MoCap model we accepted many assumptions and compromises, but gained the experience and efforts of many dedicated musculoskeletal modellers. Nonetheless, the responsibility of assessing the ability of any model to answer the question of interest rests with the researchers using the model. Our intention was to demonstrate some of the decisions and validations needed in musculoskeletal modelling of human walking and our choices reflect that.
The first of these choices is the theory on which the muscle actuation rests. The muscle models used here generate force as a function of maximum force and activation. More sophisticated muscle models, such as the Hill muscle model, are available. The Hill muscle models are thought to provide the most complete model of muscle behaviour [116], but the Hill model requires several parameters for its full definition, including: maximum isometric force, tendon slack length, optimal fibre length, pennation angle at optimal fibre length, activation and deactivation times, as well as maximum contraction velocity. The number of parameters that must be defined for a Hill muscle model, as well as the sensitivity of muscle behaviour to parameter values, means that considerable effort must be expended to parameterize and optimize Hill-type models. As with all models, poorly informed inputs (i.e. muscle parameters) lead to unstable results. This means that for certain activities it may be advisable to use a simpler, robust model rather than a more complex, unreliable one. For human activities during which muscles operate at slow contraction velocities and near their optimal fibre length (e.g. human walking), simple muscle models that depend only on geometry and a maximum force parameter are preferable [100]. Consequently, we accepted the use of the simple muscle model employed in the MoCap gait model.
Another choice was the selection of the scaling method. Understanding which scaling method to use (e.g. mass–length–fat), requires understanding the underlying assumptions of the method (e.g. distribution of mass to segments), to which model parameters the scaling applies (e.g. which segment length scale to which muscle length) and how the model uses the scaling to define parameters (e.g. are all parameters scaled or just a subset?). We made many other choices by not changing components that are inherent to the generic MoCap model, including acceptance of the form of the algorithm used to solve the muscle redundancy problem. Such choices are often embedded in standardized models, and the user does not have to modify any of them in order for the model to run to completion. This makes ensuring that the question of interest can be addressed by the model become an intentional assessment. For instance, the standard MoCap model has only one rotational DOF for the knee, which works well to assess overall walking parameters but is less well suited to a detailed study of the knee. Furthermore, standardized models are improved as more experience and data become available (e.g. [97]). For instance, the MoCap model's leg (TLEM v. 2.1) has been modified to reflect muscle paths more accurately by the inclusion of wrapping surfaces.
We also relied on the standardized definition of available outcome variables such as muscle forces and joint motions. Given that most software allows for forces and motions to be expressed in different ways (at a minimum, in different coordinate systems) it is imperative to understand the context of the reported variable. Is an output reported in the global or segment coordinate system? Is a joint force an internal force or a reaction? In what direction does the muscle force act? As with all big data, extensive data exploration and a complete understanding of the theories (mechanics), procedures (inverse dynamics) and systems (human anatomy) is necessary to draw conclusions from a musculoskeletal modelling simulation.
5.4. The role of anthropology in musculoskeletal modelling
Finally, it is important to make clear that in this demonstration we have used a walking trial of one subject and focused on the subject's left side. The kinematics and GRFs produced by this individual's right side and by other walking trials demonstrate a consistent pattern, but were not numerically identical to those reproduced here (electronic supplementary material). Consequently, the muscle forces and joint reactions are also somewhat different. The degree of intraindividual variability in gait parameters is influenced by such characteristics as age and health, but, even for this healthy adult male, peak forces within a muscle can vary. In muscles that exert relatively small forces (less than 500 N), step-to-step values can vary by more than a factor of 2 (electronic supplementary material). Muscles which show greater force generations (e.g. gastrocnemius; greater than 1000 N), step-to-step variability is lower (a factor of approx. 1–1.4; electronic supplementary material). This serves to underscore the importance of collecting and analysing multiple steps when the intent is to report patterns and forces (not the intent here).
While scaling addresses some aspects of variability and improves model predictions [35,69], scaling itself does not solve all the problems. Interindividual variability in such important variables as muscle parameters has been appreciated for over a decade [89], but the consequences of such variability for musculoskeletal modelling remains unclear. For instance, muscle strength as implemented in the solution that we use for this report depends on three factors: muscle volume, optimal fibre length and pennation angle. The MoCap model, which we used in our demonstration, uses a set of values derived from several sources but is critically informed by Klein Horsman et al. [90] and Carbone et al. [34]. These data are from measurements of cadavers who were 77 and 85 years old at the time of death, respectively. Cadaveric data are widely used in musculoskeletal models (e.g. [53,99,117]), but muscle parameters from cadavers may not represent those from individuals who are not at imminent risk of death. Table 1 includes values from Klein Horsman et al. [34] and Carbone et al. [34] and from a recent analysis of 10 healthy individuals who were evaluated with diffusion tensor imaging [91]. Gastrocnemius and soleus, the principal drivers of forward propulsion, vary in all three factors. Importantly, the complexity of gait mechanics that makes musculoskeletal modelling an attractive solution also makes understanding the effect of this variation impossible without simulating the difference. Larger volume, smaller pennation angle and shorter optimal fibre lengths all act to increase muscle strength, thereby increasing the allocation of force to it by the algorithm. But, the algorithm uses a muscle's strength relative to all other muscles' strengths, making it impossible to predict the impact of muscle parameter variation on muscle forces or joint reactions. Variation in muscle parameters is but one of the many ways that people vary in their musculoskeletal anatomy. Leveraging the expertise of biological anthropologists to assess the impact of human biological diversity on musculoskeletal modelling can only be beneficial for the entire community of researchers.
6. Conclusion
Although motivated by different specific goals, biological anthropologists and engineers (principally biomedical) both pursue questions that revolve around understanding the internal forces of the primate musculoskeletal system. Engineering questions tend to be species-specific, focusing on living humans. Biological anthropologists are also interested in living humans, but frequently ask questions that include greater taxonomic diversity. Research in both fields can be specific to the details of the individual (e.g. patient-specific orthopaedics, particular hominin fossils) or generic representatives of groups (e.g. sex-specific running shoes, comparisons of primate taxa).
Musculoskeletal modelling is powerful because of the flexibility it provides at every stage of research. Building a new model requires intensive effort, but may be necessary based on the question. Previously validated models from repositories lower barriers to modelling experiments, but researchers should be aware of choices made by model creators. Like all modelling exercises, serious and intentional thought must be dedicated to the multitude of modelling choices. This is not to say that all models need to be complex. If a simpler model adequately describes the phenomenon of interest and answers the question at hand, then a simpler model may be preferable.
Critical to support rigorous simulations, musculoskeletal models include changeable parameters that describe all aspects of the musculoskeletal system. This is a purposeful feature included by engineers that drove musculoskeletal modelling and an attractive feature for biological anthropologists. Model users can change the size and shape of the skeletal elements, dimensions and inertial properties of body segments, DOFs and range of motion of joints, action of passive soft tissue structures, the geometry and force capacity of muscles, as well as numerous other parameters. Generic models from repositories (such as the one used here) can be scaled and reshaped to capture many aspects of human variation (e.g. limb dimensions, mass), and hence are well suited to many questions. More effort can be dedicated to refine aspects of a model that are of particular interest (e.g. change the morphology of the bony pelvis). Modelling flexibility naturally extends to the information that can be derived from simulations. Forces generated by individual muscles (or portions of muscles) as well as internal to joints are standard output, as is the ability to express these quantities in different coordinate systems. Until internal forces of the body can be measured non-invasively, musculoskeletal modelling will remain the pre-eminent approach to these questions.
Future collaborations between engineers and biological anthropologists to create better models of human movement have the potential to expand the understanding of mobility in both fields. Subject-specific models are internally consistent but not generally representative of people, so extrapolating should only be done with this limitation in mind. Developing a model from scratch for a fundamentally different organism (e.g. a chimpanzee) is a significant investment.
Data accessibility
Additional information is available as the electronic supplementary material. Raw data and simulation output are available at the assigned figshare link.
Authors' contributions
A.D.S. and P.A.K. conceived the project. S.G.L. and P.A.K. designed and carried out the data collection. A.D.S. and P.A.K. carried out the simulations and data analysis. A.D.S. and P.A.K. wrote the draft and all authors contributed to the final version of the manuscript.
Competing interests
We declare we have no competing interests.
Funding
P.A.K. and S.G.L. thank S.K. Benirschke, MD, Jerome Debs Endowment Chair in Orthopaedic Traumatology at the University of Washington, for the ongoing support of this research and for our many conversations about foot form and function.
Acknowledgements
The authors would like to thank the participants of this study for their time and patience during the protocol. The authors would also like to thank the reviewers and editors for comments that improved this manuscript.
Footnotes
One contribution of 12 to a theme issue ‘Biological anthroengineering’.
Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5486386.
Review articles
A review of musculoskeletal modelling of human locomotion
,
and
Published:13 August 2021https://doi.org/10.1098/rsfs.2020.0060
Abstract
Locomotion through the environment is important because movement provides access to key resources, including food, shelter and mates. Central to many locomotion-focused questions is the need to understand internal forces, particularly muscle forces and joint reactions. Musculoskeletal modelling, which typically harnesses the power of inverse dynamics, unites experimental data that are collected on living subjects with virtual models of their morphology. The inputs required for producing good musculoskeletal models include body geometry, muscle parameters, motion variables and ground reaction forces. This methodological approach is critically informed by both biological anthropology, with its focus on variation in human form and function, and mechanical engineering, with a focus on the application of Newtonian mechanics to current problems. Here, we demonstrate the application of a musculoskeletal modelling approach to human walking using the data of a single male subject. Furthermore, we discuss the decisions required to build the model, including how to customize the musculoskeletal model, and suggest cautions that both biological anthropologists and engineers who are interested in this topic should consider.
1. Introduction
A principal research focus within biological anthropology is interpreting skeletal variation within the context of behavioural diversity, including variation in diet, disease processes and activity patterns among living and extinct primates, including humans and hominins [1–4]. Locomotor behaviour and musculoskeletal morphology have been a central focus of this area of inquiry because locomotion provides primates with access to key evolutionary resources including food, shelter and mates. Furthermore, bipedalism is considered a defining character of the hominin lineage, cementing its importance within the field [5–7]. The challenge of elucidating form–behaviour relationships are the complexities of the intervening and underlying functions (e.g. joint range of motion) and biomechanics (e.g. joint reaction forces).
Early anthropological research connecting form and function was both typological and qualitative in nature; however, anthropologists have increasingly adopted mechanical engineering approaches to answer questions about primate locomotion generally, and hominin locomotion in particular (see [4] for a historical review). Within anthropology, two major subfields of mechanics have been applied to locomotor systems: Newtonian mechanics and solid continuum mechanics. Newtonian mechanics is concerned with the description of the motion of solid bodies under the influence of a system of forces, while solid continuum mechanics quantifies the behaviour (e.g. deformation) of solid materials when subjected to forces [8,9]. Lovejoy and colleagues [10–13] were early adopters within anthropology of both Newtonian mechanics and solid continuum mechanics [4]. Newtonian mechanics is routinely used by anthropologists to investigate a variety of questions (e.g. [14–19]). Beam theory, an application of solid continuum mechanics, has been used to estimate the structural capacity of long bones and is regularly exploited by anthropologists (see [4] for a review) to reconstruct hominin locomotor adaptations (e.g. [13–17]) as well as modern human mobility patterns (e.g. [18–24]).
The application of Newtonian and continuum mechanics is, of course, a core competency of mechanical and civil structural engineering, and the methods that anthropologists use derive directly from these engineering fields. The techniques and theories that underlie their application were developed, however, to design and analyse human-made products using inert materials. While such structures and machines tend to be faithfully produced based on design specifications, humans are not clones. Rather, each individual is completed from a standard template that has, as an inherent feature, the ability to respond to the environment. This means that whereas human-made products generally have low tolerance for variance, among humans, variation is the rule. Furthermore, the ‘standard’ template itself varies among ancestral lineages, populations and species because of the evolutionary history of each (e.g. [25–27]). The biomedical engineering world has responded to this with the use of patient-specific models because, in patients, a direct relationship exists between a particular morphology (i.e. that of the patient) and the question to be answered (e.g. which implant configuration limits strain on the bone?) with the musculoskeletal model (e.g. [28–35]). Not all questions, however, concern specific individuals. In anthropology, the focus may be comparing species (or populations) to understand responses to selective pressures. In engineering, running shoes, for example, are not designed for particular individuals, who vary in foot function [36], but rather for people generally. And people are not just scaled versions of each other [37]. Thus, an appreciation of human variation, a core competency of biological anthropology, is a critical component of any modelling exercise.
At the heart of many locomotor-focused questions, both engineering and anthropological, is a need to evaluate internal (e.g. muscle, joint contact) forces. Finite-element analysis of skeletal elements [38,39], estimating locomotor costs [40,41] and evaluating task performance [42,43] all require muscle forces to be known or estimated. Measuring muscle and joint contact forces in living animals is, however, exceptionally challenging because it requires surgical intervention to implant the relevant sensors into the body [44,45]. As an additional obstacle, these techniques are normally limited to a single joint or muscle/tendon within a single study (see [44] for examples). Consequently, internal forces must be estimated. Poorly estimated loads are likely to produce misleading results, interpretations and conclusions. Consequently, estimating internal forces (e.g. muscle forces) in humans is of critical interest to both anthropologists and engineers.
1.1. Musculoskeletal modelling: a path forwards
Originally driven by a clinical interest in gait pathologies, musculoskeletal modelling has emerged as the pre-eminent technique to elucidate the internal details of human (and other animal) movement [46]. These models build upon dynamic analyses of linked rigid-segment models ([47] and described below). Rigid-segment models represent the body as a series of solid body segments which are linked together at joints. Musculoskeletal modelling extends rigid-segment models by including detailed models of individual muscles, as well as models of neurological control. In extending the dynamic approach, musculoskeletal modelling solves the muscle redundancy problem based on muscle parameters under a specified minimization criterion [48–51]. The muscle redundancy problem refers to the fact that there are more muscles than necessary to actuate movements at joints [48,49]. In solving for muscle forces, this leads to an indeterminate system of equations in which there are more unknowns (muscle forces) than constraints (equations) [51], thus requiring optimization solutions. Results of musculoskeletal modelling can reveal individual sequences of muscle excitation and activation, as well as estimates of muscle forces and bone-to-bone contact forces.
In a series of papers, Pedersen et al. (see [45]) developed an early version of this type of model and used it to estimate lower limb muscle forces and hip joint contact forces. The hip joint contact forces were validated using an instrumented hip joint replacement in a single elderly subject [45]. Consequently, these data provided an early demonstration of the validity of musculoskeletal modelling as an approach to determining internal forces. Since then musculoskeletal modelling has been used for numerous projects exploring human motion and locomotion [31,43,52–56] as well as that of non-human primates [57–59] and early hominins [60]. Wang et al. [61] applied this technique to the hominin fossil record, but their approach was necessarily limited by the lack of software available to carry out the complex calculations. Fortunately, several software packages (e.g. AnyBody Modeling System™, OpenSim, FreeBody) are now available to carry out musculoskeletal modelling [45,62,63].
Accurate musculoskeletal models demand both understanding of the principles of Newtonian mechanics and an astute appreciation of anatomy and anatomical variation. Thus, successful musculoskeletal modelling is clearly situated at the interface between engineering (i.e. mechanics) and biological anthropology (i.e. human variation). Here, we provide an overview of this effort from an anthroengineering perspective, blending knowledge from both domains with a special emphasis on a task—understanding human walking—that is of equal interest to both groups.
2. Background
2.1. The model
In musculoskeletal modelling used to solve for muscle forces, the body of interest (e.g. the whole human body) is modelled as a series of rigid segments connected by joints and actuated by muscles [34,46] (figure 1). This body moves in a virtual space whose dimensions are defined by a global Cartesian coordinate system. Models can include as few or as many body segments as necessary to answer the question of interest. For walking, models minimally include segments for the pelvis (and trunk/arms/head, often aggregated) as well as left and right feet, shanks and thighs. More complete models often also include segments that represent the patella and talus, as well as a head–arm–trunk (HAT) segment separate from the pelvic segment. Additionally, the upper part of the body may also be modelled to include segments for the lumbar region, thorax, arm, forearm, hand, neck and head segments. Body segments are defined in terms of their physical dimensions, mass, location of centre of mass and moments of inertia. Each body segment also includes its own Cartesian coordinate system, usually aligned with standard anatomical directions of the segment (e.g. proximodistal, anteroposterior, mediolateral). Other coordinate systems can be defined relative to that of the segment or to the global system and allow for ease in defining model geometry. The segments further rely on the assumption that they do not deform during locomotion (i.e. are rigid). In the software modelling/simulation environment, body segments are almost universally represented by virtual surface models of the relevant skeletal element(s).
Figure 1. MoCap model mid-simulation of walking. The image shows the model (rigid body segments represented by skeletal elements and visual representation of muscle models), four force plates (grey), the ground reaction force vector (blue line through model), force plate readings (blue lines below the force plate), experimental marker locations (blue spheres), virtual markers (red spheres with coordinate system arrows) and global coordinate system (yellow arrows).
Body segments are connected to other segments at joints, which represent articulations (e.g. hip, knee, ankle) within the body. Joints are frictionless and non-deformable and, consequently, forces and moments transfer across these joints in an equal and opposite manner. A joint definition includes its location relative to the body segments it connects, as well as the degrees of freedom (DOFs) permissible at the joint (one, two or three axes of rotations and, less commonly, translations). The direction of rotation axes are usually expressed in terms of a local (body segment) coordinate system.
With the body segments and joints defined, the skeletal portion of the musculoskeletal model is sufficient to conduct dynamic analyses. Musculoskeletal models extend the analysis by including models of individual muscles or even portions of muscles (figure 1). Muscle models range in their complexity and, as with other parts of the model, the necessary complexity depends on the motion being simulated. (For a description of how muscle complexity can be modelled, see §5.3.) A muscle model minimally requires geometry (i.e. path definition) and a parameter representing its maximum possible force. More complex muscle models also include parameters for pennation angle, optimal fibre length, physiological cross-sectional area (PCSA) and activation and deactivation times. Muscle geometry must include the coordinates of its origin (location relative to the proximal segment) and insertion (location relative to the distal segment), but muscles can be geometrically more complex than a straight line path between origin and insertion. The most basic of these complexities derives from a muscle with a broad origin (e.g. gluteus medius) or insertion (e.g. adductor magnus). In these muscles, a single origin-to-insertion path does not encompass their possible actions. One way to solve this is to split the muscle into muscle parts, referred to here as muscle elements. Each muscle element represents a region of the muscle (e.g. the superior portion of adductor magnus) and can have its own set of muscle parameters.
Most modelling systems allow for additional coordinate positions that are used to define the curved path that the physical muscles take between bony attachments (e.g. the path sartorius takes across the anterior thigh as it traverses from anterior superior iliac spine (ASIS) to proximal tibia). Via points, the simplest type of path point, are one or more coordinate positions that are intermediate between the origins and insertion points. The muscle assumes a series of straight line paths between the origin, via point(s) and insertion. The location of the via point is fixed relative to the segment's centre of mass and the connection to other segments. With via points, the muscle is required to pass through each specific point at all times. These points allow muscles to pass over other body surfaces (e.g. iliacus wrapping over the pelvic brim). In addition, muscle geometry paths can be altered by the use of ‘wrapping surfaces’ that are designed to account for relative motion between musculoskeletal elements and other soft tissue (e.g. rectus femoris sliding across the soft tissue structures anterior to the hip). These wrapping surfaces are geometric objects (e.g. cylinders, ellipsoids) that are fixed relative to a particular body segment. Muscle paths are not constrained to pass through any particular coordinate point (as with via points), but rather can move freely over the wrapping surface to find the shortest possible path across the surface between origin, insertion or via points that are beyond the wrapping surface. As the proximal and distal segments change orientation relative to each other, the point of contact between the muscle and the adjacent musculoskeletal elements can move (e.g. gastrocnemius relative to the posterodistal femur as the knee flexes). Forces developed from the action of the muscle act on the body segment used to define the geometry (origin, insertion, via points and/or points on wrapping surfaces). Muscles may also be modelled to be more complex in terms of the number of parameters that define force generation. Minimally, a muscle model must include a maximum force value which defines the force when the muscle is maximally activated, but other approximations of muscle force generation are available, e.g. the Hill-type muscle model [64].
The final required model element is a set of virtual motion markers, which is needed to match model motion to motion measured experimentally (figure 1). For motion tracking (described below), the human subject is fitted with a set of markers that are tracked by a camera system. The musculoskeletal model must include one virtual marker for each experimentally measured marker. Each virtual marker is connected to a body segment and defined in terms of that segment's coordinate system. For instance, if researchers affix a motion-tracking marker over the ASIS of a human participant, the model will require a virtual marker connected to the pelvis body segment with coordinates that position it correctly relative to the model ASIS.
2.2. The input data: motion data and ground reaction forces
An inverse dynamics approach to musculoskeletal modelling (described below) requires that the motions (translation and rotation) of the rigid segments (i.e. the skeletal elements) and externally applied forces be known. This process is well established and has been described in detail elsewhere (e.g. [65–67]). To obtain motions, retroreflective (infrared) markers attached to the skin or clothing are frequently used. Each body segment requires a minimum of three markers that are trackable throughout the entire sequence of movements, so this traceability should be checked in all portions of the data collection volume. The three-dimensional ground reaction force (GRF) for a complete stance is also necessary. GRFs are typically acquired through force plates embedded either in the floor or in a treadmill.
2.3. The simulation
2.3.1. Scaling and marker optimization
Model scaling and virtual marker position optimization are the processes of determining the model segment dimensions and virtual marker positions so that the virtual model matches the human subject and experimental protocol (motion marker placement on the human subject) as well as possible [68]. The degree of matching between the model and experimental data is evaluated using a global error metric between experimental and virtual motion marker positions. During the inverse kinematic analysis (described below), musculoskeletal modelling software finds the position of the model, including overall location in space and individual joint positions, that minimizes the sum of squared deviations between the virtual marker set and experimental motion marker data. The goal of marker optimization and body scaling is to minimize this error term.
The human models that are included with many modelling software packages are generic human models. For analysis, the model must match the subject for whom the input data are known in terms of mass, limb parameters (e.g. distances between joint centres and segment moments of inertia) and virtual marker positions [69,70]. Marker position optimization and scaling are intertwined processes because marker positions can, in part, determine segment dimensions. Initial scaling of the generic virtual model achieves a loose match between the overall bodily dimensions of the subject and those of the virtual model. If a virtual model from a repository is used (rather than creating a virtual model from scratch), then the model will have been developed using the parameters of the individual whose data were available (i.e. a generic model). This generic model will almost certainly not represent the dimensions of the subject on whom motion-tracking data and GRFs were collected, so the generic model needs to be scaled [69]. The simplest scaling approach is uniform scaling of lengths and volumes. Total body mass is the most basic dimension needed, and the ratio of the body mass of the subject to that of the generic model is typically used to adjust segment masses and other parameters that are volume-dependent. The ratio of subject to generic model stature, a length scale, is used as an initial estimation of segment dimensions and other dimensions that are length-dependent. Other anthropometric data (e.g. distances between landmarks) can be used instead of, or in addition to, stature. Other scaling methods are also available, including fully manual scaling, using body fat to adjust mass distribution among the segments, and hybrid scaling (e.g. manually scaling some parameters and algorithmically scaling others). Furthermore, depending on the sophistication of the software used to complete the mechanical analysis, these initial scalings may be adjusted during marker optimization.
Different software packages conduct marker optimization and secondary scaling of segments using different algorithms and require varying levels of interaction/input from the user. At one end of the user-input spectrum, some systems require that the user optimizes the locations of the markers on the model ‘by hand’, moving them relative to their parent segment and then interrogating error values associated with each marker during motion tracking. This process is tedious and may require significant effort before optimal marker position is achieved. At the other end of the spectrum, some software packages change both segment dimensions and virtual motion-tracking markers automatically in a single cohesive optimization process that requires little user interaction. In these frameworks, the user can permit or restrict changes in virtual marker position in none, some or all of the coordinate directions. For instance, motion-tracking markers associated with easily palpable bony landmarks (e.g. medial malleolus) are also easy to locate on the virtual model (referencing the virtual skeleton) and may be considered ‘fixed’ in their location relative to their body segment (the shank). Thus, during the optimization process, such a fixed marker would not change its position relative to the segment, but would influence the determination of shank length. Other markers (e.g. marker plates attached mid-segment) may be allowed to optimize their position in all three coordinate directions, and then contribute little to segment length determination. Finally, some markers may be well defined in two coordinate directions, but less so in the third. For instance, the superoinferior and mediolateral positions of an ASIS marker may be well defined and constrained from moving relative to the pelvis during the optimization process. The anteroposterior position of the virtual ASIS marker relative to the virtual pelvis may be more difficult to estimate as it reflects tissue thickness over this bony landmark (which is hard to measure on living subjects in the absence of internal imaging modalities). Thus, the user might fix the ASIS marker in two coordinate directions, but allow its position to optimize in the third. The fixed mediolateral and superoinferior position on the two ASISs would then influence the scaling of pelvic width.
2.3.2. Inverse kinematics/motion tracking
To carry out inverse dynamics and, in turn, estimate internal forces and moments, segment linear and rotational accelerations must be obtained. Most musculoskeletal modelling software (e.g. Anybody Modeling System™, OpenSim) uses a kinematic analysis procedure of over-determined biomechanical systems, which is formulated as a weighted least-squares optimization problem that matches experimental and model marker positions (e.g. [47,71]).
2.3.3. Inverse dynamics
While statics is the branch of mechanics that deals with stationary objects, dynamics is the branch in which motions are the focus. In statics, an equilibrium that greatly facilitates examination of the system exists where the environment (e.g. the substrate) and the structure (e.g. the human body) are balanced. In standing, the body pushes through the feet and against the ground exactly as much as the ground pushes against the body through the feet. In dynamics, the balance is time-dependent and the solution to the system is through the equations of motion, which relate forces to motions. When the forces that produce the motions are fully known, the problem is characterized as forward dynamics, while when motions are known it is described as inverse dynamics.
In human bodies, the measurement of motions (via marker tracking or other means) is much less invasive and more feasible than measuring the muscle forces that produce the motions, so currently most analyses that depend on the musculoskeletal forces generated in human gait are typically treated as inverse dynamic problems (e.g. [72–74]). Further complicating the solution in multi-body problems (i.e. those with multiple, linked rigid bodies) is that the system can be interrogated from an external or internal perspective. The external view includes the action of the GRF to produce the movement of the body along a path, while the internal one includes both the joint reactions and the action of muscle forces in creating changes in joint angles. In mechanics, the relationship between the number of unknown variables and the number of known variables determines the difficulty of the solution. Fortunately, the internal and external perspectives can be exploited to develop solutions. For instance, in the single-stance foot, the GRF is the external force and produces the moment acting on the foot segment and reacted to by the ankle joint reaction force and moment (which are an internal force and moment). Once the ankle joint reaction force and moment are known, they can be combined with the motion of the shank to solve for the knee joint reaction and moment. The inverse dynamic solution allows for this time-dependent chain solution to determine the joint reactions.
2.3.4. Muscle redundancy solution
A system where the number of variables with unknown value equals the number of independent equations that fully define the system is determinant and a unique solution exists. For systems where there are more variables with unknown values than equations—called indeterminate—assumptions are made to provide more equations [75,76]. Many biological problems, including developing muscle forces, are indeterminate because there are more muscles capable of acting across a joint than necessary from a purely mechanical perspective, i.e. some muscles are (mechanically) redundant. If the desired result from an analysis that uses inverse dynamics is the development of muscle forces, then the joint moments determined from the inverse dynamic solution need to be apportioned to individual muscle forces. This is not a trivial undertaking, because the degree of redundancy is often quite large, so no general consensus on the most accurate criterion on which to base an apportionment algorithm currently exists [77].
Musculoskeletal modelling solutions typically include an algorithm that uses muscle parameters to apportion the joint moments to particular muscles (as muscle forces). Early solutions used PCSA to solve the muscle redundancy problem [75]: the apportionment was made based on the ratio of the PCSA of a particular muscle relative to the total PCSA of all muscles capable of acting (e.g. [61]). While effective from an algorithmic perspective, PCSA-based apportionment does not reflect the natural solution, because cross-sectional area approaches assume equal activation of all possible muscles.
Another algorithm for solving the muscle redundancy problem attempts to minimize the sum of muscle activations raised to a power [76]. The activation represents a proportion of total maximum strength the muscle can generate given the current state of the model (i.e. activation is generally bound between 0 and 1). If the muscle activation exponent is 1, i.e. a linear combination, force is allocated to the muscle with the greatest moment arm until it reaches its maximum allowed value, and then other muscles are recruited. Exponents greater than 1 will distribute forces more evenly across muscles available to produce the moment at a particular joint. The greater the power, the more even will be the distribution (PCSA-based activations are essentially an exponent of infinity). Other criteria for apportionment of the joint moments have also been advanced, such as electromyography (e.g. [78,79]), but activation-based methods remain the standard [76]. Until in vivo muscle forces can be reliably measured, all solutions to the muscle redundancy problem should be treated as hypotheses about how joints and muscles produce motion.
2.4. The output
A valuable aspect of musculoskeletal modelling (or of any model) is that variables of interest can be interrogated anywhere within the model. Standard output from a musculoskeletal model includes linear and angular displacements, velocities and accelerations for all body segments. Similarly, joint positions, velocities and accelerations are also tracked and can be exported for analysis. Muscle length, activation, moment arm and force are all possible outputs, as are joint reaction forces. Additionally, users can track the position of muscle origins/via points/insertions throughout the modelled motion. Conveniently, joint reaction forces can be output in the global coordinate system or the coordinate system of the relevant body segments.
Human locomotion requires the repetition of cycles of limb and body motions that can be reduced to a stride. For walking, a step is defined from an event (e.g. initial contact) on one foot to that same event on the contralateral foot and a stride is two contiguous steps. A stride can vary in duration, even at a constant velocity, so creating an average stride or evaluating the variability of stride parameters requires normalization of the strides to some consistent duration. The typical choice is per cent of stride cycle, where the stride duration is described in terms of equal portions of the entire stride (e.g. [80,81]). Evaluating the variability of parameters of interest (e.g. muscle force, joint reactions) can be an important step in validating the calibration of a solution. While some differences among trials of the same individual are expected, an evaluation of the variability in light of the question that the model is designed to answer should occur.
2.5. Human variability
A crucial advantage of musculoskeletal modelling is the opportunity to include aspects of morphological variation, including the production of subject-specific models. Model scaling and motion tracking incorporate some aspects of human variability, including differences in body mass, limb lengths/dimensions and kinematics, but other aspects are less obviously accessible for modification in musculoskeletal modelling software. These less accessible aspects, such as variation in bone or muscle morphology, are also important to consider. Quantifying human morphological (musculoskeletal) variation is at the heart of biological anthropology. For example, human (and hominin) pelvic morphology has a profound influence on bipedal capabilities [14,80], and has been the focus of intense scrutiny [41,81,82]. Humans are known to be sexually dimorphic in pelvic size and shape [83–85]. Furthermore, human pelvic morphology is also known to vary by population [86–88]. Sources of variation can be incorporated into musculoskeletal modelling to understand pathology [31] and differences in joint reaction forces [35]. For example, marker-informed parameter optimization and secondary scaling (described above) can change the width of the pelvis if experimental and virtual markers are included on bony landmarks of the pelvis, such as the ASISs. If the subject is wider than the generic model, the parameterization can widen the pelvis so that the virtual and experimental markers are more closely aligned. Scaling can also be accomplished manually either by applying factors to the three segment coordinate directions or by redefining the locations of specific anatomical landmarks associated with the segment. The former method moves all points proportionately while the latter allows for localized warping.
In addition to skeletal variation, aspects of muscle morphology that vary among humans may also be critical for accurate musculoskeletal modelling. For instance, variability in muscle parameters has been of concern for some time [89]. While repository models represent a particular person's morphology, people differ in important musculoskeletal modelling parameters such as muscle volume, optimal fibre length and pennation angle (table 1).
Table 1. Comparison of muscle parameters between the standard MoCap model in AMMR™ v. 2.3 (AnyBody Technology, Aalborg, Denmark; from Klein Horsman et al. [90] and Carbone et al. [34]) and Charles et al. [91]. Muscle names are per the MoCap AMMR v. 2.3. Where multiple muscle regions are present in the MoCap model, the values for the muscle volumes are summed while the optimal fibre lengths and pennation angles are averaged.Collapsemuscleoptimal fibre length (cm)pennation angle (°)muscle volume (ml)MoCapCharles—averageMoCapCharles—averageMoCapCharles—average
adductor brevis | 10.38 | 7.6 | 0.0 | 11 | 108.90 | 93 |
adductor longus | 10.57 | 11.0 | 0.0 | 12 | 159.56 | 159 |
adductor magnus | 9.97 | 23.1 | 0.0 | 12 | 559.25 | 567 |
biceps femoris caput breve | 9.14 | 10.9 | 0.0 | 9 | 107.95 | 92 |
biceps femoris caput longum | 8.54 | 20.4 | 29.9 | 11 | 232.01 | 194 |
extensor digitorum longus | 5.99 | 13.8 | 8.3 | 7 | 32.29 | 76 |
extensor hallucis longus | 5.98 | 10.6 | 14.4 | 7 | 36.27 | 21 |
flexor digitorum longus | 3.84 | 28.4 | 25.28 | |||
flexor hallucis longus | 5.20 | 30.1 | 161.50 | |||
gastrocnemius lateralis | 5.69 | 12.2 | 25.4 | 9 | 136.36 | 128 |
gastrocnemius medialis | 6.01 | 9.7 | 10.8 | 10 | 263.26 | 230 |
gemellus | 3.43 | 0.0 | 28.40 | |||
gluteus maximus | 13.57 | 0.0 | 936.55 | |||
gluteus medius | 6.76 | 5.3 | 674.71 | |||
gluteus minimus | 3.28 | 0.0 | 82.59 | |||
gracilis | 18.11 | 17.3 | 0.0 | 6 | 87.97 | 91 |
iliacus | 8.16 | 8.8 | 203.13 | |||
obturator externus inferior | 6.88 | 0.0 | 37.88 | |||
obturator externus superior | 2.77 | 0.0 | 68.18 | |||
obturator internus | 2.05 | 0.0 | 52.08 | |||
pectineus | 9.50 | 0.0 | 64.58 | |||
peroneus brevis | 4.54 | 23.1 | 86.09 | |||
peroneus longus | 5.08 | 15.8 | 121.52 | |||
peroneus tertius | 4.27 | 19.1 | 26.52 | |||
piriformis | 3.88 | 0.0 | 31.25 | |||
plantaris | 4.77 | 0.0 | 11.36 | |||
popliteus | 2.40 | 7.4 | 0.0 | 8 | 25.57 | 15 |
psoas major | 9.92 | 13.4 | 193.18 | |||
psoas minor | 5.91 | 0.0 | 6.63 | |||
quadratus femoris | 3.37 | 0.0 | 49.24 | |||
rectus femoris | 7.84 | 14.2 | 21.9 | 8 | 226.33 | 249 |
sartorius | 34.71 | 40.8 | 0.0 | 205.49 | 143 | |
semimembranosus | 8.09 | 15.8 | 25.0 | 12 | 138.26 | 244 |
semitendinosus | 14.16 | 18.3 | 0.0 | 8 | 208.33 | 186 |
soleus | 4.40 | 14.6 | 61.6 | 12 | 792.91 | 461 |
tensor fasciae latae | 9.49 | 0.0 | 83.33 | |||
tibialis anterior | 6.83 | 13.7 | 9.6 | 7 | 181.49 | 129 |
tibialis posterior | 3.78 | 34.2 | 163.36 | |||
vastus intermedius | 7.68 | 18.1 | 11.8 | 11 | 292.61 | 521 |
vastus lateralis | 6.69 | 19.6 | 0.0 | 15 | 583.33 | 606 |
vastus medialis | 7.16 | 15.9 | 0.0 | 14 | 454.27 | 415 |
Several possibilities for future work that explores the impact of morphological variation on muscle and joint forces are worth noting. For instance, while some anthropological work has elucidated the interactions of morphology and walking conditions (e.g. velocity, changes of direction, burdens and gradients) on energy expenditure (e.g. [92–94]), much less is known about how morphology and walking condition factors interact to affect joint and muscle forces (although see [31] for an example). People vary morphologically and walking under varied conditions is an activity of daily living, so this area is an important one for many disciplines to pursue including biological anthropology, biomedical engineering and product design. Perhaps more importantly, collecting the requisite data to inform a musculoskeletal model is time consuming, invasive and dependent on (at this time) fragile laboratory equipment. Consequently, acquiring samples that represent the breadth of human morphological variation is unlikely to occur in the near future. Once validated, though, a musculoskeletal model can be modified to assess the impact of these known morphological differences on muscle and joint forces. Extending that line of research further, a longer-term goal would be use musculoskeletal modelling to understand the influence of the morphological differences between modern humans and extinct hominins on muscle and joint forces, an area that has sought biomechanical solutions without effective software support for decades [12,14,61,95,96].
2.6. Current study
Here, we demonstrate the process and output of one human subject during walking using AnyBody Technology's musculoskeletal modelling software (Anybody Modeling System™; AnyBody Technology, Aalborg, Denmark). We present several of the possible outputs that could be of interest to biological anthropologists. This includes joint kinematics, joint reaction forces and a selection of lower limb muscle force magnitudes. We make some recommendations for future researchers to help facilitate their use of software and output. Additionally, by presenting a set of possible outputs, we demonstrate the motivation for undertaking musculoskeletal modelling and that results may have a direct bearing on anthropological questions, or as intermediate steps to other analyses (e.g. finite-element analysis). Additionally, while engineering has much to offer biological anthropologists, the latter have an appreciation for human variation which may enhance future generations of musculoskeletal models.
3. Material and methods
3.1. Participant
One healthy participant was used for this demonstration, but he is part of a larger study aimed at examining human walking. The participant was 36 years old (male; body mass, 74.3 kg; stature, 172 cm) and reported being free from lower-limb injuries. The University of Washington's Institutional Review Board approved all aspects of this study (IRB no. STUDY00001125).
3.2. Experimental protocol
Kinetic and kinematic data were measured using a four-force plate (Kistler, Switzerland), 10-camera motion capture system (Qualisys, Sweden) in the Amplifying Movement & Performance (AMP) laboratory of the University of Washington. Thirty infrared-retroreflective 5 mm hemispherical markers were affixed to anatomical landmarks and used to track motion through the laboratory (figure 1; electronic supplementary material). The participant walked unshod the length (10 m) of the gait laboratory five times at his self-selected preferred walking speed (1.15 m s−1). The participant walked in a straight path with his eyes directed at a point on the far wall of the laboratory and contacted the surface of each force plate with a single foot. Trials with multiple feet on a single force plate or ones where the stance foot crossed plates were immediately redone. The participant was allowed several attempts prior to data collection to become familiar with the protocol. All trials exhibited similar kinematics and kinetics (data not shown), so we selected one trial to use in this demonstration.
Marker data and GRFs were collected at 120 Hz and 1200 Hz, respectively, and then were filtered at 5 Hz and 10 Hz, respectively, with a fourth-order, low-pass zero-lag Butterworth filter. Calibration of the system yielded a limitation in its fidelity for marker data of 1 mm and force data of ±2.5 N for the direction of travel (X), ±5 N side (Y) and ±25 N vertical (Z). All data were exported from the Qualisys Track Manager software in C3D file format, which can be read directly by the modelling software (provided in electronic supplementary material).
3.3. Musculoskeletal model
We used the ‘MoCap’ model from the AnyBody Managed Model Repository™ (AMMR v. 2.3; AnyBody Technology, Aalborg, Denmark; http://www.anybodytech.com/software/model-repository-ammr), which is a multi-trial, full-body, motion-capture-driven human gait model, and the commercially available AnyBody Modelling System™ (v. 7.3; AnyBody Technology, Aalborg, Denmark) to calculate forces in the lower limb [34,97]. The MoCap model that we used for this demonstration has been used to assess human gait in numerous applications (e.g. [98,99]). The MoCap model is assembled from a trunk and left and right lower limb components. The trunk includes a lumbar region, trunk, neck and head. The lower limbs consist of the following segments: thigh, patella, shank, talus and foot. Each lower limb has six total DOFs, including all three rotations at the hip and one each at the knee (flexion/extension), ankle (plantar flexion/dorsiflexion) and subtalar (inversion/eversion) joint. The pelvis (relative to the ground) has six DOFs, three translational and three rotational. Forty-one anatomically defined, bilateral muscles (table 1) are represented by 169 muscle elements in each model lower limb (e.g. gluteus medius is composed of 12 separate muscle element actuators) [34]. The muscle elements in this model generate force as a function of their activation and a strength parameter which is appropriate for human walking [97,100]. We defined a virtual motion-tracking marker for each of the experimentally measured motion markers (electronic supplementary material). The polynomial algorithm that was employed to solve the muscle redundancy problem apportioned force to muscles using an objective function with power = 3. (For more information regarding the model, see §5.3.)
3.4. Musculoskeletal simulation
The first step of the simulation was to scale the model to match the dimensions of the subject while simultaneously optimizing the marker coordinate and joint rotation axes' locations [68]. After initial scaling, the values of these parameters (i.e. segment dimensions, marker locations, rotation axes' locations) were determined using the marker motion data from a trial. Values for all parameters were optimized to minimize the sum of square distances between model marker positions and experimental marker positions. Following convergence of parameter optimization, we conducted an inverse kinematics analysis. Finally, we conducted the inverse dynamic analysis, which includes apportioning muscle forces using the algorithm that minimizes the cost function as a sum of muscle activation values raised to the third power. This value has been shown to allocate muscle activations in a way that reflects experimental data [77,100,101].
3.5. Output variables and data processing
We output several kinematic and kinetic variables, including joint kinematics, filtered GRF components and joint reaction forces (electronic supplementary material). Additionally, we output the force magnitude for a selection of the muscle elements for the left lower limb during an entire stride cycle of walking. Muscle element forces within anatomically defined muscles were summed for presentation (e.g. 12 elements that represent gluteus medius). All time-related curves were resampled to 1% increments of the stride cycle.
4. Results
GRF component curves across the gait cycle are presented in figure 2. As is typical for human walking, the vertical component has two peaks (early stance peak, 719 N; late stance peak, 778 N) on either side of a valley (633 N) which occurs around the middle of stance phase. The average of the two peaks is 103% of body weight (729 N), while the valley represents 87% of body weight.
Figure 2. Filtered ground reaction force components. Grey vertical line represents moment of toe-off.
Pelvic rotation and obliquity and hip, knee and ankle kinematics are provided in figure 3. All kinematics are of the left lower limb throughout a single stride. Pelvic motion is described in terms of the right hip (previous stance limb) relative to the left hip (current stance limb) [102]. Hip, knee and ankle joint reaction forces, which include contributions from muscle forces, are provided in figure 4. In each joint, the vertical joint reaction force is greater in magnitude than the anteroposterior and mediolateral components. All three joints show a large peak during the second half of the stance phase, and this late peak is significantly larger in the ankle. Muscle force profiles for 18 anatomically defined muscles across the stride cycle are provided in figure 5. During the stance phase of walking, gastrocnemius, soleus and gluteus medius generate the largest forces. All output data are provided in the electronic supplementary material.
Figure 3. Angular kinematics for the pelvic segment and hip, knee and ankle joints. Vertical grey line represents toe-off.
Figure 4. Joint reaction forces expressed in the global coordinate system. Vertical grey line represents toe-off.
Figure 5. Muscle forces for 18 anatomically defined muscles.
5. Discussion
5.1. Model output
The subject that we used to demonstrate the musculoskeletal modelling approach to modelling walking is an adult male without gait pathology and within species-typical anthropometric characteristics. Consequently, the model-derived motions should accord with those determined experimentally and extensively described in the last 40 years of gait research (e.g. [66,103–105]). The model results presented above do meet this requirement, as discussed below in detail. We focus on the model-derived kinematics, because these kinematics are the foundation for all subsequent analyses in musculoskeletal models. Any analysis should begin by verifying that model-derived variables align with those obtained empirically, and the cause of any discrepancies should be understood or noted as a potential limitation. We belabour the details in the following paragraphs in order to demonstrate the basic approach to model validation.
In our subject at heel strike, the pelvis is rotated approximately 10° in the horizontal plane such that the right hip is posterior to the left hip (e.g. the current stance limb). The pelvis then swings through a neutral position and, while leading up to the next right heel strike, into an orientation with the right hip anterior to left hip (approx. 10°). The pelvis drops in the frontal plane with the right hip inferior to the left (stance) hip soon after heel strike as the right lower limb moves into the swing phase with pelvic obliquity ranging between ± 4°. The pelvis remains in this position until the subsequent right heel strike, after which the pelvis approaches the neutral position. Model-derived pelvic rotation and obliquity are similar to empirically determined values (e.g. [103]).
In the sagittal plane, the left hip begins the stance phase in flexion of approximately 17° and moves through a neutral position and into extension of approximately 17°. Winter [66] describes a less symmetrical flexion–extension cycle at natural cadence (mean values of approx. 20° of flexion and 11° of extension), but the model-derived hip rotations are well within one standard deviation of his mean values. During most of the stance phase, the left thigh is adducted less than 6° in accord with others (e.g. [103]).
While the pelvic and hip joints have three rotational DOFs, the knee and the ankle only rotate in the sagittal plane. The left knee flexes slightly (10°) early in the stance phase, but maintains a near fully extended position through midstance, after which it begins to flex in preparation for the swing phase where flexion exceeds 60°. After the initial heel strike and as the forefoot contacts the substrate, the ankle plantar flexes slightly (less than 5°) and then assumes a dorsiflexed position for most of stance. Rapid plantar flexion occurs in late stance as the body is propelled forwards, reaching a plantar flexion value of greater than 10°. The ankle returns to a dorsiflexed position for swing. These values for both knee and ankle angles are typical (e.g. [66,103,105]).
Consequently, the kinematics of the model match those that we expect empirically. While this check is not sufficient to ensure that the internal (muscle and joint reaction forces) are reasonable, it is a necessary condition and validation for any model. Incorrect kinematics guarantee incorrect or suspect muscle forces.
Joint reaction forces have been empirically measured only with instrumented joint implants. By definition, the surgery creates conditions that limit applicability of the results, but no other method currently exists to measure force. Pedersen et al. [45] measured hip joint contact forces in a 72-year-old man (63.2 kg) with a hip replacement and found that the vertical force was 1750 N, considerably less than the almost 3500 N we found with the musculoskeletal model. Important differences between this patient and our subject exist, including age (72 versus 36 years), body mass (62.3 versus 74.3 kg), health status (58 days post hip replacement surgery versus healthy) and walking speed (0.89 versus 1.15 m s−1). While hip and knee joint reaction forces presented here are larger than in vivo values using instrumented prostheses [77,106,107], our model-derived joint forces are, however, similar to other models of healthy adults [35,107].
Although it is difficult to compare muscle forces with experimental data, it is possible to consider the pattern of activation with reference to the muscle's action. The muscles that generate the highest force in the model are gluteus medius, soleus and gastrocnemius. Gluteus medius plays a critical role in bipedal walking by preventing pelvic drop during single stance [103,105]. The gluteus medius muscle profile from the model is double peaked, and the timing of the two peaks is coincident with the two peaks in the vertical GRF, indicating the muscle's support of the pelvis during single stance. Gastrocnemius and soleus plantar flex the ankle; when the foot is in contact with the ground and the hip is extended, plantar flexing the ankle moves the body anteriorly [103,105]. The gastrocnemius and soleus muscles show the greatest force generation during the propulsive portion (second half) of the stance phase, in keeping with our expectation from their function.
Other muscles also display moderately high forces. Tibialis anterior, which dorsiflexes the ankle, shows muscle force generation immediately after heel strike when the body begins to move over the foot, and then is mostly inactive until the second half of the swing phase. The quadriceps demonstrate a force peak early in stance phase, and again at the stance–swing transition, potentially relating to their role in extending and stabilizing the knee [105]. The hamstring muscles show the greatest peak at the end of swing phase and into early stance phase. Activation of these muscles at the end of swing phase is associated with decelerating the limb prior to heel strike. Early in stance the GRF is directed anterior to the hip and knee, and the hamstring muscles are active to prevent further flexion (hip) and extension (knee) of these joints. The hip flexors are active late in the stance phase, and may be active to counteract hamstring hip extension moment as the hamstrings flex the knee in preparation for swing phase. The adductors show initial peak in stance (adductor magnus) and late stance (adductor longus), but generally smaller than forces in other muscles.
5.2. Input data choices
The model-derived motions and joint reaction and muscle forces align with our expectations from the experimental data in the literature. We achieved these results through careful attention to detail when we collected the data we used as inputs to the musculoskeletal model. We provide below some brief suggestions from our experience. Because motions are derived from the change in position of these markers, marker placement on the segment should be chosen to define fully the motion of the segment and to reduce marker positional noise as much as possible. Positional noise arises from the inherent error of measurement of the marker's position in space (system error) and from movement of the surface (skin or clothing), where the markers are placed relative to the underlying rigid segment (e.g. [108–110]). System error is reduced through calibration; for our system, it is less than 1 mm for any marker. Movement between the external (skin) surface and the bone is more complicated to reduce, but may not be clinically relevant [110].
Marker position is a compromise with soft tissue coverage of underlying bone, landmark palpability (to aid placement) and distance between markers (such that marker positional error is much less than the distance between markers) all being important. The shank and foot have palpable landmarks that are covered by minimal soft tissue (e.g. the malleoli), but the thigh (e.g. greater trochanter) and, especially, the pelvis are more challenging segments on which to affix markers. This problem is exacerbated when the subjects of interest have more soft tissue (muscle or fat) overlying the bony landmark, and considerable soft tissue coverage is common in the hip–pelvis region. One approach to ameliorate some of this difficulty is to use a marker plate which has four (or more) markers mounted on a hard plastic plate that is then strapped to the segment. If a plate can be affixed firmly to the segment, as it can be to the thigh, the marker plate eliminates between-marker error. Unfortunately, marker plates are difficult to affix firmly to the pelvis. Marker placement issues, consequently, are most difficult to overcome to track the motion of the pelvis. Yet, for musculoskeletal models of human locomotion, pelvic motions are critical components of system performance. We achieved the motions shown in figure 3 through extensive protocol development and testing.
Another critical input to the simulation are the GRFs. It is possible to conduct an inverse dynamics solution of the single support phase of walking with a single force plate (e.g. [111]), but extending the analysis to the double support portions of the stride is more difficult with one force plate because it would require estimation of the unknown contribution of the other stance foot. Consequently, any research that requires information about double support should optimally be conducted with multiple force plates. The placement of the force plates relative to each other should be appropriate to the population and gait. Stride lengths for adults or running gaits are longer than for children or walking, and plates should be positioned to prevent subjects ‘targeting’ foot placement or large numbers of repeated trials. Data analysis is less complicated if stance occurs on one plate, but the data from two immediately adjacent plates can be used as long as there is no gap between them [112]. Multiple feet contacting a single plate is also possible to resolve with additional information (such as localized pressure sensors) or predictive methods (e.g. [113,114]). All of these remediation techniques are time-intensive. The laboratory in which the presented data were collected has four floor-embedded force plates that are positioned to allow normal stride lengths for adults when walking. The subject and trial that we used for this demonstration is part of a much larger study [115]. Nonetheless, we assessed each foot placement on all four plates for each trial before moving on to the next trial. Such data checking is a tedious effort at the time of data collection, but well worth the effort to ensure all data are suitable for future analyses. Using four force plates, we had data frames before and after the data that we used in the analysis. This allowed us to discard the initiation and termination frames of the simulation, which are frequently unusable. Consequently, researchers are advised to locate force plates well and test the placement prior to data collection.
Finally, researchers may consider collecting complete anthropometrics of subjects and using these in the initial scaling as the starting points for the parameterization. People vary considerably in terms of their proportions [35,37], so standard models, such as the model we used here, represent a morphology that may not match the dimensions of any particular subject. If your model starts with segment lengths that are substantially different from the segment lengths of the model, the parameterization algorithm may not be able to converge on a suitable solution. Two failures can occur: the parameterization can fail to converge, which produces an obvious error message, or the parameterization converges, but, in order to converge, the simulation requires changes to the morphology or motions that are inconsistent with the experimental data or subject dimensions. This latter failure can be insidious. An absolutely critical step in model validation is to check the motions and shapes of the segments after parameterization. A poorly parameterized model is incorrect and will not produce acceptable motion or joint and muscle forces.
5.3. Modelling choices
In order to produce this demonstration, we made a number of choices with regard to the model which bear discussion. The most important of these was the choice to use a standard, repository model, in this case the MoCap model of AMMR v. 2.3. By using the MoCap model we accepted many assumptions and compromises, but gained the experience and efforts of many dedicated musculoskeletal modellers. Nonetheless, the responsibility of assessing the ability of any model to answer the question of interest rests with the researchers using the model. Our intention was to demonstrate some of the decisions and validations needed in musculoskeletal modelling of human walking and our choices reflect that.
The first of these choices is the theory on which the muscle actuation rests. The muscle models used here generate force as a function of maximum force and activation. More sophisticated muscle models, such as the Hill muscle model, are available. The Hill muscle models are thought to provide the most complete model of muscle behaviour [116], but the Hill model requires several parameters for its full definition, including: maximum isometric force, tendon slack length, optimal fibre length, pennation angle at optimal fibre length, activation and deactivation times, as well as maximum contraction velocity. The number of parameters that must be defined for a Hill muscle model, as well as the sensitivity of muscle behaviour to parameter values, means that considerable effort must be expended to parameterize and optimize Hill-type models. As with all models, poorly informed inputs (i.e. muscle parameters) lead to unstable results. This means that for certain activities it may be advisable to use a simpler, robust model rather than a more complex, unreliable one. For human activities during which muscles operate at slow contraction velocities and near their optimal fibre length (e.g. human walking), simple muscle models that depend only on geometry and a maximum force parameter are preferable [100]. Consequently, we accepted the use of the simple muscle model employed in the MoCap gait model.
Another choice was the selection of the scaling method. Understanding which scaling method to use (e.g. mass–length–fat), requires understanding the underlying assumptions of the method (e.g. distribution of mass to segments), to which model parameters the scaling applies (e.g. which segment length scale to which muscle length) and how the model uses the scaling to define parameters (e.g. are all parameters scaled or just a subset?). We made many other choices by not changing components that are inherent to the generic MoCap model, including acceptance of the form of the algorithm used to solve the muscle redundancy problem. Such choices are often embedded in standardized models, and the user does not have to modify any of them in order for the model to run to completion. This makes ensuring that the question of interest can be addressed by the model become an intentional assessment. For instance, the standard MoCap model has only one rotational DOF for the knee, which works well to assess overall walking parameters but is less well suited to a detailed study of the knee. Furthermore, standardized models are improved as more experience and data become available (e.g. [97]). For instance, the MoCap model's leg (TLEM v. 2.1) has been modified to reflect muscle paths more accurately by the inclusion of wrapping surfaces.
We also relied on the standardized definition of available outcome variables such as muscle forces and joint motions. Given that most software allows for forces and motions to be expressed in different ways (at a minimum, in different coordinate systems) it is imperative to understand the context of the reported variable. Is an output reported in the global or segment coordinate system? Is a joint force an internal force or a reaction? In what direction does the muscle force act? As with all big data, extensive data exploration and a complete understanding of the theories (mechanics), procedures (inverse dynamics) and systems (human anatomy) is necessary to draw conclusions from a musculoskeletal modelling simulation.
5.4. The role of anthropology in musculoskeletal modelling
Finally, it is important to make clear that in this demonstration we have used a walking trial of one subject and focused on the subject's left side. The kinematics and GRFs produced by this individual's right side and by other walking trials demonstrate a consistent pattern, but were not numerically identical to those reproduced here (electronic supplementary material). Consequently, the muscle forces and joint reactions are also somewhat different. The degree of intraindividual variability in gait parameters is influenced by such characteristics as age and health, but, even for this healthy adult male, peak forces within a muscle can vary. In muscles that exert relatively small forces (less than 500 N), step-to-step values can vary by more than a factor of 2 (electronic supplementary material). Muscles which show greater force generations (e.g. gastrocnemius; greater than 1000 N), step-to-step variability is lower (a factor of approx. 1–1.4; electronic supplementary material). This serves to underscore the importance of collecting and analysing multiple steps when the intent is to report patterns and forces (not the intent here).
While scaling addresses some aspects of variability and improves model predictions [35,69], scaling itself does not solve all the problems. Interindividual variability in such important variables as muscle parameters has been appreciated for over a decade [89], but the consequences of such variability for musculoskeletal modelling remains unclear. For instance, muscle strength as implemented in the solution that we use for this report depends on three factors: muscle volume, optimal fibre length and pennation angle. The MoCap model, which we used in our demonstration, uses a set of values derived from several sources but is critically informed by Klein Horsman et al. [90] and Carbone et al. [34]. These data are from measurements of cadavers who were 77 and 85 years old at the time of death, respectively. Cadaveric data are widely used in musculoskeletal models (e.g. [53,99,117]), but muscle parameters from cadavers may not represent those from individuals who are not at imminent risk of death. Table 1 includes values from Klein Horsman et al. [34] and Carbone et al. [34] and from a recent analysis of 10 healthy individuals who were evaluated with diffusion tensor imaging [91]. Gastrocnemius and soleus, the principal drivers of forward propulsion, vary in all three factors. Importantly, the complexity of gait mechanics that makes musculoskeletal modelling an attractive solution also makes understanding the effect of this variation impossible without simulating the difference. Larger volume, smaller pennation angle and shorter optimal fibre lengths all act to increase muscle strength, thereby increasing the allocation of force to it by the algorithm. But, the algorithm uses a muscle's strength relative to all other muscles' strengths, making it impossible to predict the impact of muscle parameter variation on muscle forces or joint reactions. Variation in muscle parameters is but one of the many ways that people vary in their musculoskeletal anatomy. Leveraging the expertise of biological anthropologists to assess the impact of human biological diversity on musculoskeletal modelling can only be beneficial for the entire community of researchers.
6. Conclusion
Although motivated by different specific goals, biological anthropologists and engineers (principally biomedical) both pursue questions that revolve around understanding the internal forces of the primate musculoskeletal system. Engineering questions tend to be species-specific, focusing on living humans. Biological anthropologists are also interested in living humans, but frequently ask questions that include greater taxonomic diversity. Research in both fields can be specific to the details of the individual (e.g. patient-specific orthopaedics, particular hominin fossils) or generic representatives of groups (e.g. sex-specific running shoes, comparisons of primate taxa).
Musculoskeletal modelling is powerful because of the flexibility it provides at every stage of research. Building a new model requires intensive effort, but may be necessary based on the question. Previously validated models from repositories lower barriers to modelling experiments, but researchers should be aware of choices made by model creators. Like all modelling exercises, serious and intentional thought must be dedicated to the multitude of modelling choices. This is not to say that all models need to be complex. If a simpler model adequately describes the phenomenon of interest and answers the question at hand, then a simpler model may be preferable.
Critical to support rigorous simulations, musculoskeletal models include changeable parameters that describe all aspects of the musculoskeletal system. This is a purposeful feature included by engineers that drove musculoskeletal modelling and an attractive feature for biological anthropologists. Model users can change the size and shape of the skeletal elements, dimensions and inertial properties of body segments, DOFs and range of motion of joints, action of passive soft tissue structures, the geometry and force capacity of muscles, as well as numerous other parameters. Generic models from repositories (such as the one used here) can be scaled and reshaped to capture many aspects of human variation (e.g. limb dimensions, mass), and hence are well suited to many questions. More effort can be dedicated to refine aspects of a model that are of particular interest (e.g. change the morphology of the bony pelvis). Modelling flexibility naturally extends to the information that can be derived from simulations. Forces generated by individual muscles (or portions of muscles) as well as internal to joints are standard output, as is the ability to express these quantities in different coordinate systems. Until internal forces of the body can be measured non-invasively, musculoskeletal modelling will remain the pre-eminent approach to these questions.
Future collaborations between engineers and biological anthropologists to create better models of human movement have the potential to expand the understanding of mobility in both fields. Subject-specific models are internally consistent but not generally representative of people, so extrapolating should only be done with this limitation in mind. Developing a model from scratch for a fundamentally different organism (e.g. a chimpanzee) is a significant investment.
Data accessibility
Additional information is available as the electronic supplementary material. Raw data and simulation output are available at the assigned figshare link.
Authors' contributions
A.D.S. and P.A.K. conceived the project. S.G.L. and P.A.K. designed and carried out the data collection. A.D.S. and P.A.K. carried out the simulations and data analysis. A.D.S. and P.A.K. wrote the draft and all authors contributed to the final version of the manuscript.
Competing interests
We declare we have no competing interests.
Funding
P.A.K. and S.G.L. thank S.K. Benirschke, MD, Jerome Debs Endowment Chair in Orthopaedic Traumatology at the University of Washington, for the ongoing support of this research and for our many conversations about foot form and function.
Acknowledgements
The authors would like to thank the participants of this study for their time and patience during the protocol. The authors would also like to thank the reviewers and editors for comments that improved this manuscript.
Footnotes
One contribution of 12 to a theme issue ‘Biological anthroengineering’.
Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5486386.