The development of a segment-based musculoskeletal model of the lower limb: introducing FreeBody.

Traditional approaches to the biomechanical analysis of movement are joint-based; that is the mechanics of the body are described in terms of the forces and moments acting at the joints, and that muscular forces are considered to create moments about the joints. We have recently shown that segment-based approaches, where the mechanics of the body are described by considering the effect of the muscle, ligament and joint contact forces on the segments themselves, can also prove insightful. We have also previously described a simultaneous, optimization-based, musculoskeletal model of the lower limb. However, this prior model incorporates both joint- and segment-based assumptions. The purpose of this study was therefore to develop an entirely segment-based model of the lower limb and to compare its performance to our previous work. The segment-based model was used to estimate the muscle forces found during vertical jumping, which were in turn compared with the muscular activations that have been found in vertical jumping, by using a Geers' metric to quantify the magnitude and phase errors. The segment-based model was shown to have a similar ability to estimate muscle forces as a model based upon our previous work. In the future, we will evaluate the ability of the segment-based model to be used to provide results with clinical relevance, and compare its performance to joint-based approaches. The segment-based model described in this article is publicly available as a GUI-based Matlab® application and in the original source code (at www.msksoftware.org.uk).

Introduction

Musculoskeletal modelling technology promises to provide an estimation of the forces experienced by the hard and soft tissues of the body during movement [1-4]. This is important given both the difficulties associated with in vivo measurements of these forces and the potential for improved treatment outcomes if the result of an intervention could be simulated a priori.

One approach to musculoskeletal modelling of the lower limb is to predict the muscular, ligamentous and joint reaction forces that create movement from measurements of the movement outcome (i.e. the kinematics of the body segments and the ground reaction force (GRF)), an approach which is referred to as inverse dynamics. Traditional approaches to inverse dynamics tend to pose the equations of motion in terms of the forces and moments that act across the joints. This simplifies the description of motion by making it unnecessary to explicitly describe all the forces and moments experienced by the limb which, instead, are implicitly captured by the concept of the joint. The alternative to this ‘joint-based’ approach is to pose the equations of motion in terms of the forces and moments that act upon the body segments (a ‘segment-based’ approach) which requires all the forces and moments to be explicitly described. We have recently demonstrated that segment-based approaches can provide additional insights as to the functional anatomy of the lower limb [5,6], not least because of the additional detail that is incorporated within the approach. However, these analyses have been based on the consideration of simple two-dimensional models and not a detailed three-dimensional (3D) model of the lower limb.

To date, few 3D segment-based models have been proposed within the literature (notable exceptions being the lower limb model of Moissenet et al. [7] and the upper limb model of Pennestri et al. [8]). In our previous work, we have described a detailed, 3D, optimization-based approach to inverse dynamics analysis which we have argued has advantages over traditional approaches [9-11] (these findings have been supported by other groups [7,12]). However, this model is not entirely segment-based in its approach. In particular, a number of joint-based assumptions are made in order to simplify the modelling of the knee. The purpose of this study was therefore to develop a fully segment-based model of the lower limb and to compare its predictions to our previous work. Finally, the model described within this paper has been made publicly available, and this is also described towards the end of this paper.

Material and methods

This study describes the development of a new musculoskeletal model of the lower limb that is based upon our previous work [9-11,13]. The lower limb is modelled as a 3D arrangement of five rigid segments that represent the foot, shank, thigh, pelvis and patella (the patella segment is assumed to have zero mass). The data inputted into the model consists of motion capture data describing the position of retro-reflective markers attached to a subject (kinematic data) and ground reaction forces (kinetic data). The motion capture data are used to specify the position and orientation of the segments as a function of time, and from this the kinematics of the segments can be calculated. The model is then used to estimate the muscular and joint reaction forces experienced by the lower limb during the recorded movement. The forces that create the observed movement are therefore predicted from their resultant movements and forces (inverse dynamics analysis).

Model posture

The locations and orientations of the foot, shank, thigh and pelvis segments are derived from raw data that describe the location (in the laboratory fixed global coordinate frame; GCS) of retro-reflective markers placed upon key anatomical landmarks (table 1). Each segment has a segment fixed local coordinate frame (LCS) which is defined by reference to the position of the markers when the subject is stood in a known calibration position (the anatomical position; figure 1). Within each LCS, the y-axis is first defined to run from the distal joint to the proximal joint. An intermediate z-axis from medial to lateral is then defined by reference to anatomical landmarks, and then the x-axis found such that it is orthogonal to the y- and intermediate z-axes. Finally, the permanent z-axis is found such that it is orthogonal to the x- and y-axes. The origin of each LCS is located at the centre of the proximal joint (apart from the pelvis, whose origin is at the midpoint of a line connecting the anterior superior iliac spines). In constructing the LCS of each segment, the transformation of the segment from the LCS to the calibration position is found. Next, the transformation of each segment from the calibration position to its position in each frame is found using the method of Horn [14]. Ultimately, the location and orientation of a segment for each frame is defined by the transformation from the LCS of the segment to its position in the GCS (this can be found by combining the former two transformations). Within the GCS, the x-axis runs from posterior to anterior, the y-axis from inferior to superior and the z-axis from medial to lateral.

Table 1.

Marker positions used for data capture.

marker	location
FCC	calcaneus
FMT	tuberosity of the fifth metatarsal
FM2	head of the second metatarsal
TF	additional marker placed on the foot
FAM	apex of the lateral malleolus
TAM	apex of the medial malleolus
C1, C2, C3	additional markers placed on the shank segment
FLE	lateral femoral epicondyle
FME	medial femoral epicondyle
T1, T2, T3	additional markers placed on the thigh segment
RASIS	right anterior superior iliac spine
LASIS	left anterior superior iliac spine
RPSIS	right posterior superior iliac spine
LPSIS	left posterior superior iliac spine

Figure 1.

Position of markers on a subject standing in the calibration (anatomical) position.

Patellar posture

The position and orientation of the patella is given as a function of the tibiofemoral joint flexion angle (θ) using a new model of the patella that has been developed from the literature as follows. First, the location of the patella origin is established relative to the position of the tibial tuberosity in the tibial LCS (this is taken from the data of Klein Horsman et al. [15] and scaled in the same way as is described later for the other muscle parameters). This is achieved by determining the path of the patellar tendon and by defining the origin of the patella to be at the point where the patellar tendon meets the patella. Specifically, the orientation of the patellar tendon relative to the tibial LCS can be described by the combination of the patellar tendon sagittal and coronal plane angles. These are calculated from the knee flexion angle using polynomial regression equations derived from the results of Kobayashi et al. [16] (table 2). Combining these angles with the length of the patellar tendon (which is taken from Klein Horsman et al. and scaled to the subject characteristics) allows the position of the patella origin to be found.

Table 2.

Coefficients describing the position and orientation of the patella as a function of knee flexion angle. (The applicable equation is of the form: variable= b0+b1θ+b2θ2+b3θ3+b4θ4, where θ is the knee flexion angle. PT, patellar tendon.)

variable	b0	b1	b2	b3	b4
PT sagittal plane angle	20.4	−2.60×10−1
PT coronal plane angle	10.9	−2.33×10−1	1.89×10−3	−5.69×10−6
patellar flexion	5.59	6.60×10−1
patellar tilt	1.63	6.67×10−2	1.44×10−4	−5.37×10−6
patellar rotation	1.43	1.06×10−1	−3.45×10−3	5.47×10−5	−2.38×10−7

Next, the orientation of the patella is established, again by reference to the knee flexion angle. The orientation of the patella is given based upon a set of three Euler angles. The sequence of rotation is as follows: (i) patellar flexion about the z-axis (of the femoral LCS); (ii) patellar tilt about the new y′ axis; and (iii) patellar rotation about the new x′′ axis. These three angles are also calculated using polynomial regression equations, this time developed from the data of Nha et al. [17]

The above calculations allow the transformation from the patellar LCS to the GCS to be established. This in turn means that the position and orientation of each segment for each time point has been established. Note that anatomical landmarks were used to combine the above datasets in order to ensure their compatibility.

Adding the musculoligamentous geometry

Once the locations and orientations of every segment for every frame have been determined, it is then possible to add the next layer of detail—that is the geometry of the muscular and ligamentous structures. This is achieved by using the data of Klein Horsman et al. [15] which specify the origin, insertion and path of 163 muscle elements and 14 ligaments. The data of Klein Horsman et al. is used to define the locations of the origins and insertions of the muscles and ligaments within the LCS of the applicable segments (this is achieved by constructing an LCS for each of the Klein Horsman et al. segments using exactly the same methodology as is used to define the LCS of the subject). The geometry of the musculoligamentous system can then be calculated for each frame using the transformations from LCS to GCS described earlier.

Scaling

The musculoskeletal geometry of Klein Horsman et al. [15] is scaled to match the size of the subject under consideration. This is achieved at the segmental level, by calculating the relative size of the subject's segments in comparison to that of the Klein Horsman cadaver, and then adjusting the coordinates of origins and insertions accordingly (linearly within the Cartesian 3D frame).

Wrapping

For muscle and ligament elements that do not wrap around other musculoskeletal structures, their line of action is taken to be the straight line from origin to insertion. For those elements that do wrap around other structures, the wrapping of the element is represented in one of two ways. For the majority of the muscles, the wrapping is described by the description of additional ‘via’ points through which the path of the muscle is constrained to pass, and that thus define the line of action of the muscle. These via points are also taken from the work of Klein Horsman et al. However, in the case where the muscle element is free to move over the underlying surface, a wrapping cylinder is defined to represent the underlying structures, and then the path of the element around the cylinder is found using the method of Charlton & Johnson [18]. The wrapping cylinders are also taken from the work of Klein Horsman et al., however, an additional cylinder was also defined to represent the wrapping of the quadriceps tendon around the femoral condyles at deeper knee flexion angles [19].

Anthropometry and kinematics

The anthropometry of the model (i.e. the mass, centre of mass and inertia tensor of each segment) is generated using the data of de Leva [20]. The kinematics of the movement (i.e. the linear and angular velocities and accelerations of each segment) is calculated using the quaternion-based methodology of Dumas et al. [21].

Calculating the model kinetics

The analysis approach used in this study follows the example of our previous work [9-11]. That is, the indeterminate equations of motion governing the movement of the lower limb (in the GCS) are posed based upon the musculoskeletal geometry and the measured kinematics (segment motions) and kinetics (ground reaction forces) of the lower limb. An optimization approach is then used to solve the equations of motion simultaneously. In this study, a number of different formulations of the equations of motion are employed (table 3), in order to evaluate the effect of transitioning to a segment-based model and adding more detail relating to the geometry of the knee joint, and these are described below. It should be noted that in all cases the method of Dumas et al. [21] is used to formulate the equations of motion.

Table 3.

Details of the different cases considered in this study.

case	patellofemoral joint modelled?	two tibiofemoral joint contacts?	joint reaction forces explicit?
1	no	no	no
2	yes	no	no
3	yes	no	yes
4	yes	yes	yes

Case 1 comprises the same solution approach as presented in our previous work [11], although using the revised musculoskeletal geometry described above. First, the inter-segmental forces can be found directly by the iterative application of Newton's second law to each segment, working from distal to proximal: Next, the equations of motion are posed by considering the motion of foot, shank and thigh segments. This allows 18 indeterminate equations of motion to be written (equation (2.2); table 4)—nine describing the linear motion of the segments (three for each segment) and nine describing the rotational motion (again three for each segment). The equations of motion are parametrized by the previously calculated musculoskeletal geometry, the kinematics of the segments and the GRF. The 186 unknowns consist of the muscular, ligamentous and joint reaction forces (bone on bone contact forces) experienced by the lower limb. It should be noted that each segment has a full complement of 6 degrees of freedom (d.f.) and there are no constraints imposed based upon considerations of the kinematics of the joints (with the exception of the patellofemoral joint). This means that each joint effectively also has 6 d.f. (apart from the patellofemoral joint), although the joints are not explicitly described in the model: In equation (2.2), the rotational effect of each muscle and ligament element upon a segment k is captured by calculating an effective moment arm ( and , respectively) which describe the overall rotational effect of 1 N of tension in the element on the given segment [5]. So for a monoarticular muscle, includes the effect of both the muscle force and the joint reaction force that arises because of the muscle force. This is akin to the assumptions made in a joint-based approach. Equally, for the intermediate segment of a biarticular muscle that does not attach to the segment, captures the rotational effects of the joint reaction forces created by tension in the muscle [10].

Table 4.

Nomenclature used in equations of motion.

a^^k	linear acceleration of the centre of mass of segment k
c^^k	vector from centre of rotation of joint at proximal end of segment k to centre of mass of segment k
d^^k	vector from centre of rotation of joint at proximal end of segment k to centre of rotation joint at distal end of segment k
d~^k	skew-symmetric matrix of vector
d~l3	skew-symmetric matrix of vector from centre of rotation of hip to tibiofemoral joint contact l
E3×3	3×3 matrix of zeros
f~^3	skew-symmetric matrix of vector from centre of rotation of hip to contact point of patella with the femur
Fi	magnitude of force in muscle i
F_max_i	maximum possible force in muscle i (upper bound)
g^	acceleration due to gravity
h~l2	skew-symmetric matrix of vector from centre of rotation of knee to tibiofemoral joint contact l
i	muscle number
I3×3	3×3 identity matrix
j	ligament number
J	cost function
k	segment number
Lj	magnitude of force in ligament j
L_maxj	maximum possible force in ligament j (upper bound)
mk	mass of segment k
M	total number of muscles
N	total number of ligaments
p^ik	unit vector representing the line of action of force created by muscle i that acts on segment k (zero if muscle does not insert on segment k)
pat	patella
pt	patellar tendon
P/Q ratio	ratio of patellar tendon to quadriceps tendon force
q^jk	unit vector representing the line of action of force created by ligament j that acts on segment k (zero if ligament does not insert on segment k)
r^ik	vector from centre of rotation of joint at proximal end of segment k to point of action of muscle i on segment k (zero if muscle does not insert on segment k)
Rxk	x component of reaction force acting at proximal end of segment k
Ryk	y component of reaction force acting at proximal end of segment k
Rzk	z component of reaction force acting at proximal end of segment k
R^^k	vector representing x, y and z components of reaction force acting at proximal end of segment k
R^lk	vector representing x, y and z components of reaction force l acting at proximal end of segment k
s^jk	vector from centre of rotation of joint at proximal end of segment k to point of action of ligament j on segment k (zero if ligament does not insert on segment k)
−S^^k	inter-segmental force acting on proximal end of segment k
v^ik	effective moment arm of muscle i on segment k
w^jk	effective moment arm of ligament j on segment k
−W^^k	inter-segmental moment acting on proximal end of segment k
Y3×3k	inertia tensor of segment k
ρi	P/Q ratio for muscle i (zero if the muscle is not part of the quadriceps muscle group)
φ^˙^k	angular velocity of segment k
φ^¨^k	angular acceleration of segment k

In case 1, the effect of tension in the patellar tendon is assumed to create an equal and opposite effect on the tibial and femoral segments, and the patellofemoral joint is not explicitly modelled (this is again a common assumption in a joint-based approach). However, the patella is assumed to be in force and moment equilibrium and thus the P/Q ratio (the ratio of the patellar to quadriceps tendon forces) will alter with knee flexion [22]. The P/Q ratio is derived from a consideration of the musculoskeletal geometry of the patella in the sagittal plane (figure 2) and then the effective force upper bound of the patellar tendon is adjusted as a function of the knee flexion angle to reflect this relationship.

Figure 2.

Moment equilibrium of the patella results in a changing ratio between the patellar and quadriceps tendon forces (P and Q, respectively) which depends on the angles of incidence of the two tendons on the patella (p and q, respectively).

In case 2, an explicit consideration of the patella is used to derive the equations of motion. First, there is the addition of three more unknowns representing the x, y and z components of the patellofemoral joint reaction force. The patellofemoral joint reaction force is included in the equations of motion of the thigh segment, and in addition, the lines of action of the quadriceps muscle elements are explicitly modelled. Second, three more equations of motion are added that describe the force equilibrium of the patella. Finally, a fourth equation of motion is added which ensures that the sagittal plane components of the patellar tendon and quadriceps forces are in the same ratio as that which was derived from geometric considerations in case 1. It should be noted that the patella is also assumed to be in moment equilibrium but this assumption is not explicitly captured in the equations of motion. Instead, the pivot point of the patella is assumed to shift in order to maintain equilibrium. The patellofemoral joint contact position is found for each frame by finding this pivot point, which is in turn used to define the contact position between the patella and thigh segments.

The equations of motion in case 2 (equation (2.3); table 4) thus consist of 21 equations and 190 unknowns (three additional unknowns representing the patellofemoral joint reaction force and one additional unknown representing the patellar tendon force which is now modelled separately from the quadriceps components): For cases 1 and 2, the joint reaction forces (apart from the patellofemoral joint reaction force in case 2) are not described explicitly in the moment parts of the equations of motion. The approach is therefore a hybrid between joint- and segment-based approaches. In case 3, the joint reaction forces are also explicitly included in the moment parts of the equation of motion (equation (2.4); table 4). Otherwise, case 3 is the same as case 2. The position at which the joint reaction forces are considered to act is defined in the following way. The joint contact points are assumed to be the joint centres given in the Klein Horsman et al. [15] dataset. These are defined to be fixed within the distal segment of a joint—thus the contact position on the proximal segment is not fixed, and will change if there is a translation of the two segments relative to one another. The contact surfaces are assumed to be rigid in this model:

In case 4, the tibiofemoral joint contact force is compartmentalized into a medial and a lateral component (equation (2.5); table 4). The number of unknowns representing the tibiofemoral joint contact force is thus increased from 3 (x, y and z components) to 6. Thus, the number of unknowns for case 4 is 193.

For cases 1–3 the joint reaction forces are considered to act through the centres of rotation of the joints taken from the Klein Horsman et al. dataset [15]. In case 4, the lateral and medial tibiofemoral contact forces act through points 2 cm laterally and medially to the centre of rotation (an inter-condyle distance of 4 cm was chosen following the example of previous work [23]). The tibiofemoral joint contact force created by each muscle that spans the knee must be split between these lateral and medial compartments. This cannot be achieved by a consideration of the force parts of the equation of motion alone. This provides a further reason why it is important to explicitly include the joint reaction forces within the moment parts of the equations of motion: For each case, the equations of motion were solved for each time point using an optimization approach that is based upon a cost function adapted from the work of Crowninshield & Brand [24] and Raikova [25] and that we have used previously [11] (equation (2.6); table 4). The optimization was solved using the optimization toolbox of Matlab® (R2103a; The Mathworks, Inc, 2013):

The force upper bound of each muscle () was calculated from the data of Klein Horsman et al. [15]. Specifically, the cross-sectional area of each muscle element was doubled to account for the larger muscle mass of our young, athletic population and then multiplied by an assumed maximum muscle stress (3.139×105 N m−2). [26] The force upper bounds of the ligaments () were similar to our previous work [11] (table 5). It should be noted that the ligaments are simply modelled as tensile force actuators. That is, the calculated ligament forces are not derived from measurements of ligament strain. A detailed justification and analysis of this approach (including a discussion of limitations) is provided in our previous work [11].

Table 5.

Upper bounds of the ligaments included in this study.

ligament	joint	upper bound (N)
iliofemoral ligament (anterior)	hip	850
iliofemoral ligament (lateral)	hip	850
pubofemoral ligament	hip	450
ischiofemoral ligament	hip	450
anterior cruciate ligament	knee	2000
posterior cruciate ligament	knee	4000
medial collateral ligament	knee	3000
lateral collateral ligament	knee	2000
oblique popliteal ligament	knee	1000
posterior tibiotalar ligament	ankle	850
tibiocalcaneal ligament	ankle	850
tibionavicular ligament	ankle	850
posterior talofibular ligament	ankle	850
calcaneofibular ligament	ankle	850

Experimental data

The specific data used in this study has also been previously described [13,27,28] and is the same data that was used in our previous work describing optimization approaches to inverse dynamics analysis [9-11] (this being considered an advantage as it makes comparison between studies much easier). To summarize, the data describe the performance of vertical jumps performed by a group of athletic males (n=12; mean age 27.1±4.3 years; mean mass 83.7±9.9 kg) who provided informed consent prior to the data collection process (the data collection was approved by the institutional review board of St Mary's University College). The position of retro-reflective markers placed upon key anatomical landmarks [29,30] was captured at 200 Hz using Vicon (Vicon MX System, Vicon Motion Systems Ltd, Oxford, UK) and ground reaction forces were recorded synchronously using a Kistler force plate (Kistler Type 9286AA, Kistler Instrumente AG, Winterthur, Switzerland). Both the position and force data were filtered using a fifth order Woltring filter [31] before it was input into the model. The filtered raw data are provided with this article as the electronic supplementary material.

Statistical analysis

The performance of each case was evaluated by comparing the predicted muscle forces to electromyography (EMG) envelopes, representing the muscular activations during vertical jumping. These envelopes were taken from the work of Rodacki et al. [32] (who generated these envelopes from the raw electromyograms by using the Myo-Dat v. 5.0 EMG analysis software (MIE Medical Research Ltd., Leeds, UK) with a second-order low pass filter of 6 Hz). First, the mean muscle force at each time point (relative to the time of take-off) was calculated in order to produce a composite curve of muscle forces for each case. The similarity of this curve to the EMG data were quantified by using a Geers' metric (equation (2.7)) [33], which it has been suggested is appropriate for comparing measured and experimental curves in biomechanics [34-36]. The Geers' metric provides separate estimations of the magnitude error (M, which is insensitive to differences in phase between the two curves) and the phase error (P, which is insensitive to differences in magnitude between the two curves): where and m(t) and c(t) are the measured and calculated waveforms, respectively.

Results

On average, the optimization was able to find a solution for 95.1±6.2% of the analysed frames. The lowest overall joint contact forces were found in case 4 (table 6). There was a marked similarity in the joint contact forces found in cases 2 and 3, and the ankle and total tibiofemoral joint contact forces found in case 1 were also similar to cases 2 and 3.

Table 6.

Mean peak joint contact forces (mean ± s.d.) predicted for each case in this study (TFJ, tibiofemoral joint contact force; PFJ, patellofemoral joint contact force).

		knee
case	ankle	lateral TFJ	medial TFJ	total TFJ	PFJ	hip
1	6.7±1.1			8.5±1.6		8.7±2.0
2	6.7±1.1			8.5±1.6	9.9±3.0	7.2±1.9
3	6.7±1.1			8.5±1.6	10.0±3.0	7.3±1.8
4	6.7±1.1	2.7±1.3	4.5±0.6	6.2±1.0	8.2±2.5	6.2±1.3
Cleather et al.[11]	9.0±1.6			7.4±2.1		5.5±1.0

There was a marked similarity between the predicted muscle forces and the EMG data across cases, as indicated by the similarity in the mean Geers' metrics (tables 7 and 8). For instance, figure 3 illustrates the mean muscle forces of six major muscle groups for case 4, in comparison to the EMG envelopes of Rodacki et al. [32]. There was a clear qualitative similarity between the model predictions and the EMG envelopes. The Geers' metric suggested a close agreement for soleus, rectus femoris, glutaeus and vastus, but demonstrated a lower level of agreement for gastrocnemius and the biarticular hamstrings.

Table 7.

Magnitude error (Geers' metric) for the comparison of average predicted muscle forces with EMG envelopes for each case considered in this study (gastroc., gastrocnemius; r. fem., rectus femoris; bi. ham., biarticular hamstrings).

case	soleus	gastroc.	r. fem.	vastus	glutaeus	bi. ham.	mean ± s.d.
1	−0.05	−0.51	−0.11	−0.22	−0.19	−0.27	−0.22±0.14
2	−0.05	−0.49	−0.02	−0.25	−0.24	−0.30	−0.23±0.16
3	−0.08	−0.53	−0.08	−0.25	−0.24	−0.30	−0.25±0.15
4	−0.07	−0.50	−0.11	−0.25	−0.22	−0.33	−0.25±0.14

Table 8.

Phase error (Geers' metric) for the comparison of average predicted muscle forces with EMG envelopes for each case considered in this study (gastroc., gastrocnemius; r. fem., rectus femoris; bi. ham., biarticular hamstrings).

case	soleus	gastroc.	r. fem.	vastus	glutaeus	bi. ham.	mean ± s.d.
1	0.11	0.33	0.07	0.14	0.13	0.21	0.16±0.09
2	0.11	0.34	0.04	0.15	0.15	0.22	0.17±0.09
3	0.12	0.34	0.05	0.16	0.15	0.23	0.17±0.09
4	0.12	0.32	0.07	0.15	0.15	0.25	0.18±0.08

Figure 3.

Average muscular recruitment during vertical jumping as predicted by case 4. The solid line represents the average percentage of maximum force capability expressed during jumping. The shaded grey regions represent EMG envelopes that have been reported previously in the literature (Rodacki et al. [32]). M is the magnitude error and P is the phase error (both taken from the Geers' metric).

Discussion

This study describes the development of a new musculoskeletal model of the lower limb that has been developed from our previous work in the estimation of muscle and joint contact forces during vertical jumping. The motivation for this was to create a musculoskeletal model that is segment-based. The results of this study demonstrate that the segment-based model has a similar ability to estimate muscle forces during vertical jumping as the previous approaches that we had described.

The performance of the model was assessed by comparing the muscle force estimations to muscular activations that had been previously reported within the literature. The mean muscle force estimations for each case were compared with the mean activation (EMG envelopes) using a Geers' metric, which provides values representing both the magnitude and phase errors. The differences between the estimated forces and the EMG envelopes showed only minor variations among cases, and the models demonstrated the closest match for soleus and rectus femoris, and the biggest discrepancies for gastrocnemius and the biarticular hamstrings (tables 7 and 8). Taken as a whole, the analysis presented here provides new and quantitative evidence as to the validity of this lower limb model.

The Geers' metric (phase error) demonstrated that the estimated muscle forces showed a consistent time lag relative to the EMG envelopes. Some degree of time lag would be expected because of the neuromechanical delay between muscular activation and force production, and so a zero phase error would not represent the best fit, and the phase error here is in the expected direction. Equally, the Geers' metric (magnitude error) demonstrated that the magnitudes of the model estimates were consistently lower than the muscular activations. A major factor in this trend was that the duration of peaks in estimated muscle forces was lower than the duration of the activations. It should be noted that muscle estimations and EMG envelopes that were compared in this study were not recorded concurrently, nor taken from the same subjects. Although this is a clear limitation of the study, and could account for some of the differences between the curves, there is evidence that the muscular activations during vertical jumping are remarkably stable between subjects and across jump conditions [32]. Therefore, given that the comparisons here are between group mean curves this was considered an acceptable limitation. In future work, we will seek to make the same comparisons at the subject level.

It is notable that the estimated biarticular muscle forces (biarticular hamstrings and gastrocnemius) showed a lower level of agreement than the monoarticular muscle force estimates. This was predominantly because of a later onset of muscular activation in the estimated muscle forces than was seen in the EMG envelopes. We have previously identified that modelling the function of the biarticular muscles is a key challenge for the musculoskeletal modelling community [9,10]. In particular, the benefits derived from the biarticular muscles transferring power between joints [37] may not be entirely captured by this modelling approach.

The modelling approach described in this paper is markedly different to other contemporary models of the lower limb. Many of these differences are consistent with, permitted by, or a consequence of the segment-based approach to biomechanics that is used here. Some of the more important differences and their ramifications are described below.

One trend within musculoskeletal modelling is to optimize the measured kinematics based upon kinematic constraints that describe the joints [7,12,38] (although there are other kinematic optimizations that do not assume joint constraints). One advantage of this is that these constraints can be used to limit the effects of measurement errors (for instance, soft tissue artefacts). However, the optimization in turn reduces the d.f. of the model. This may be disadvantageous, as these d.f. could be important in describing the mechanics of the limb. For instance, it is common to pose lower limb models which do not permit joint translations, however, we have recently demonstrated that the anterior–posterior translation of the tibiofemoral joint surfaces contributes to the ability of the musculature of the knee to rotate the thigh [5]. Equally, we have argued that it is more appropriate to retain the full amount of d.f. [39] and to instead represent the joints by modelling the forces that constrain the kinematic behaviour of the joint [1]. We believe that this is important as some of the structures that provide the forces that constrain the kinematics are those that are commonly injured.

The model that is described in this article does not impose such kinematic constraints. Instead, the position and orientation of each segment is determined independently. The potential problem of soft tissue artefacts is ameliorated by using a redundant number of markers on each segment, and then using the method of Horn [14] to find the best fit. All d.f. are consequently captured within this model, such that joint translations are permitted and all joints have three rotational d.f. The larger number of d.f. does bring its own costs however. In particular, the model presented here is not able to find a solution for all of the frames analysed. We believe this is a problem that would be eliminated as more detail is added to the model (by providing a greater number and variability in the lines of action of force to compensate for the greater number of d.f. [1,39-41]).

A common simplification in musculoskeletal modelling is the employment of multiple step solution processes that consider the equations of motion related to forces and moments independently [1,7,11]. The advantage of these approaches is that they reduce the complexity of the equations of motion, making their solution more straightforward. However, such an approach is not physiologically realistic—in the body forces and moments are equilibrated simultaneously and this dynamic interplay is important in understanding how the tissues of the body are loaded. In our previous work, we have described a one step, simultaneous inverse optimization approach which permits a more faithful representation of the function of biarticular muscles [9,10] and musculoligamentous interaction [11] during movement, and the same approach is employed in this model. Latterly, some other groups have adopted a similar approach and found it to be advantageous [7,12]. In particular, Moissenet et al. [7] have shown how such an approach can be used to successfully estimate the timing of musculotendon forces during gait. Equally, Hu et al. [12] found that such an approach predicted joint contact forces during gait that more closely matched the forces measured in patient populations than less complex (more traditional) approaches. Preliminary work therefore seems to suggest that more detailed modelling of the interplay of forces in creating linear and rotational motion is important to the understanding of movement.

One example of a contemporary multi-step approach is in the partitioning of tibiofemoral joint contact forces into medial and lateral components. For instance, Gerus et al. [23] describe a process whereby first the joint moments are calculated, second the muscle forces are determined based upon the calculated joint moments and then finally the joint contact forces are found such that the sum of the moments of muscle and joint contact forces is equal to the joint moment. This process is repeated twice, once for each of the lateral and medial compartments. However, this methodology has the same weaknesses as the multi-step approaches previously referred to. In particular, the medial and lateral tibiofemoral contact forces contribute to the equilibration of moments at the knee, and thus should be solved for simultaneously with the muscle and ligament forces. Only in this approach is the dynamic interplay that equilibrates forces and moments properly represented. In the model described here, the tibiofemoral joint contact forces are found simultaneously with the muscle and ligament forces (this represents a considerable advance over our previous work [9-11]).

When employing a joint-based, multi-step approach it is common to not include an explicit representation of the patellofemoral joint. This is an understandable limitation of the approach—if the lower limb is modelled as a chain of linked segments, and analysed in terms of the inter-segmental moments between those segments, it is difficult to include a consideration of the patella. However, in order to create a segment-based description of the lower limb it was necessary for us to develop a patella model and to include the patellofemoral joint contact force. The patella model incorporated within this model is simple—the movement of the patella is entirely determined by the flexion of the tibiofemoral joint, and this relationship is not changed to reflect subject-specific characteristics. The movement of the patella can have an important impact on the outputs of a musculoskeletal model, in large part due to its effect on the P/Q ratio and in turn the tension required in the quadriceps musculature. In our future work, we will develop the model described here to permit a subject-specific scaling of the patellofemoral model which, in particular, produces an accurate, subject-specific P/Q ratio.

To summarize, this paper has described the development of a segment-based model of the lower limb and has demonstrated that its muscle force estimates are in line with our previous work. The rationale for such a model was based on our recent theoretical work that demonstrated advantages of a segment-based approach [5,6]. However, this study has not sought to evaluate whether the segment-based model presented here has such advantages in comparison to more traditional joint-based models. One reason for this is that such a comparison is very difficult to make—there are a number of assumptions that are variable between the two approaches. Despite this, we intend to address this question in our future work, and the development of this model is an important step along the way.

Our future plans are not limited to the segment-/joint-based comparison. Work is already underway to try and validate the muscle force estimates of the model by comparison with directly measured experimental data and values from the literature [35]. We have also performed some preliminary work evaluating the ability of the model to detect clinically relevant differences in tibiofemoral joint contact forces following an exercise intervention [42]. In addition, we believe that the lack of understanding as to the behaviour of musculoskeletal models is a key impediment to the realization of a clinical tool [1,43]. It is our intention to perform a systematic, probabilistic analysis of the input parameters and assumptions behind FreeBody. This will be invaluable in guiding researchers towards those parameters which are most likely to require accurate subject-specific information, and the degree to which subject-specific accuracy is required. In addition, the statistical analysis will provide a robust estimate as to the potential error inherent to the current state of the art approaches.

Introducing FreeBody

Currently, the choice of musculoskeletal modelling software available for use is somewhat limited. The market is dominated by a small number of key products which include AnyBody (AnyBody Technology A/S, Aalborg, Denmark), SIMM (MusculoGraphics Inc., Santa Rosa, CA, USA), LifeModeler (LifeModeler Inc., San Clemente, CA, USA) and OpenSim (Simbios, Stanford, CA, USA). There is no doubt that these software packages represent the state of the art in musculoskeletal modelling software, especially given that they tend to be backed by some of the most illustrious research groups within the field (e.g. SIMM and OpenSim are associated with the research group of Scott Delp at Stanford University and AnyBody with John Rasmussen of Aalborg University). However, for the most part these packages are commercial concerns, and carry large licence fees which could be prohibitive for some users. One exception is the freely available OpenSim [44], which is used quite widely within the biomechanics community.

A further barrier to the wider adoption of musculoskeletal modelling software is that such packages tend to require a considerable time investment to learn about their operation and then to develop applications with the desired functionality. In particular, this will serve as a barrier to potential users who do not have a strong background in biomechanics, software programming, mathematics or physics. There is therefore a need for a musculoskeletal modelling application that is more user friendly and that makes this analysis approach available to a wider population of potential users.

The musculoskeletal model described in this article as case 4 is publicly available as both a Matlab application and in the original source code (at www.msksoftware.org.uk) and includes extensive documentation. In addition, the source code for the model and the Matlab application are provided with this article as the electronic supplementary material. The Matlab application is driven by a graphical user interface which makes the use of the model straightforward and intuitive and the ubiquity of Matlab within science and education means that many users will be very familiar with the model environment. We, therefore, hope that this version is almost entirely ‘plug and play’, that is that the user simply has to process their data into the appropriate input format, and then they can use the model as an analysis tool. It is our intention that this brings musculoskeletal modelling technology to an entirely new population of users. Our philosophy is also to create a research tool that has broad use for the widest possible range of users, and it is for this reason that we have released the underlying code. Ultimately, it is our hope that this release of FreeBody will help to encourage the development of technology within this area.