Estimation of musculoskeletal models from in situ measurements of muscle action in the rat hindlimb

Although studies of biological motor control often focus on the strategies implemented by neural systems, the properties of peripheral systems – the skeleton, muscles, tendons, connective tissue, etc. – also make significant contributions to the production of movement. In many cases these peripheral systems have been shown to simplify neural control or make natural behaviors more efficient (Nishikawa et al., 2007; Pfeifer et al., 2007; Valero-Cuevas et al., 2007b). Any understanding of the neural control of movement, therefore, clearly requires an understanding of the peripheral musculoskeletal systems with which it interacts.

Such an understanding of peripheral systems is commonly formalized in computational musculoskeletal models (e.g. Delp and Loan, 1995; Sueda et al., 2008; Teran et al., 2003; Valero-Cuevas, 2009). Musculoskeletal models have been used in many different contexts and have been helpful in gaining insights into biological motor control. Although there are a wide variety of musculoskeletal models used in the literature, in general these models are created from ‘the bottom up’: investigators make a large number of independent measurements of physiological parameters (e.g. optimal muscle length, physiological cross-sectional area, moment arms, etc.) which are then used to specify a model of muscle action (Brown and Loeb, 2000; Burkholder and Nichols, 2004; Holzbaur et al., 2005; Johnson et al., 2008; Kargo et al., 2002; Kargo and Rome, 2002; McKay et al., 2007; Topp and Rosen, 1990). Such reductionist approaches are very appealing because they provide a standardized way to construct musculoskeletal models across a wide range of species, and measured parameters are usually directly related to well-defined muscle properties.

Although this approach to developing musculoskeletal models has been generally successful, there are several issues concerning this approach that are often difficult to evaluate. At a basic level, the validity of such models can be difficult to establish: the measurement of muscle parameters is subject to error and this error can result in a range of estimated muscle actions (Santos et al., 2009; Santos and Valero-Cuevas, 2003; Valero-Cuevas et al., 2003). One crucial property of muscles included in such models is their moment arm. In muscles that are represented as line segments, the moment arm of a muscle depends on its origin and insertion. However, many muscles have a broad insertion and/or origin. If a muscle is represented with a single origin and insertion to define its line of action, it is not always obvious, based on anatomical dissection, where to best place these points. Furthermore, it is not always clear what degree of model complexity is necessary to capture muscle action. For instance, to model especially complex muscles with very broad insertions or origins, investigators often divide the muscle into several distinct regions, modeling each region with an independent insertion and origin which together span the observed muscle anatomy (Chadwick et al., 2008; Holzbaur et al., 2005). How to divide these complex muscles into distinct regions and the most appropriate location of their respective insertion and origin points is not usually obvious.

The present study develops and evaluates an approach that can potentially help guide and validate the development of musculoskeletal models. We take an integrative approach to modeling the actions of individual muscles, using ‘top down’ empirical characterizations of measured muscle actions to create parameterized muscle models. This approach is an extension of that developed by Loeb et al. (Loeb et al., 2000) and is related to other work attempting to estimate muscle moment arms by direct measurement of muscle action (Lee et al., 2008; Pearlman et al., 2004; Valero-Cuevas et al., 2000). In particular, we measured, in situ, the forces and moments produced by stimulation of individual muscles in the rat hindlimb. These forces and moments are a direct consequence of the intrinsic force-generating properties of a muscle working through its attachment to the skeleton (Loeb et al., 2000). Any muscle model attempting to characterize the action of a muscle should therefore be able to predict these forces and moments well. In addition, because these measurements are closely linked to the properties of the stimulated muscle, it should be possible to estimate the parameters for a muscle model (i.e. its anatomical origin and insertion and its intrinsic force-generating properties) based on the forces and moments observed experimentally. The present study took this approach, using the experimentally observed forces and moments to create a parameterized muscle model to capture the action of muscles.

We applied this approach to two different muscles in the rat hindlimb, the gracilis posticus and the biceps femoris (Greene, 1935). The gracilis posticus is a relatively simple, band-like muscle with well-defined origin and insertion. We expected that this muscle would be well described with a model consisting of a single insertion and origin. In contrast, the biceps femoris is a complex muscle with an insertion distributed across the tibia from the knee to the ankle (Gillis and Biewener, 2001). Given this complex insertion, it is not obvious where one should place the insertion of this muscle on the tibia, or whether multiple insertions might be required to describe the action of this muscle. We expected that this muscle would be poorly described by a single insertion and origin and would require a more complex model structure to capture its action (Valero-Cuevas et al., 2007a). These two muscles, representing a range from simple to complex muscle architectures allowed us to evaluate the utility of our approach for the development and validation of muscle models.

A preliminary report of these experiments and analyses was previously published as a conference proceeding (Tresch et al., 2010).

Data was obtained from five adult female rats (Ratus norvegicus Berkenhaut 1769; mass 250–300 g). All procedures were approved by the Institutional Animal Care and Use Committee at Northwestern University. Rats were anesthetized with ketamine and xylazine (80 and 10 g kg^–1; intraperitoneally). Metal posts were attached to three points on the pelvis using bone screws and cement: one each, on the rostral and caudal ends of the contralateral (right) pelvis and another on the rostral end of the ipsilateral pelvis (Fig. 1A). A threaded attachment was then cemented to bone screws attached at the distal tibia so that it was approximately perpendicular to the tibia. The anterior head of gracilis was severed and gently peeled back to expose the nerve innervating the posterior head of gracilis. The nerve was gently freed from the surrounding tissue and a small cuff electrode with bipolar exposed leads was placed around the nerve, insulated with Vaseline, and secured to surrounding tissue using cyanoacrylic glue. To prepare the biceps femoris (BF) for stimulation, the anterior head of the biceps femoris (BFa) was first severed close to its origin on the vertebrae and partially removed by gently cutting its adhesion to the posterior portion of the biceps (BFp). BFp here refers to the entire region of the BF, which has an origin on the ischium and inserts along the extent of the tibia (see Discussion). We then exposed the branch of the tibial nerve innervating the BFp by cutting the gluteus maximus and caudofemoralis, then peeling back the caudofemoralis to expose the nerve. The branch of the nerve to the semimembranosus was cut; the remaining nerve was freed from surrounding tissue, placed in a cuff electrode, insulated with Vaseline, and secured to the surrounding tissue.

We then separated the semitendinosus from the BF, cutting through the connections between the two muscles and severing the semitendinosus distal to the fusion of its medial and lateral heads. Exposed muscles were covered with moist gauze. Finally, we performed a dorsal laminectomy to expose the lumbar spinal cord and dorsal roots. The L4 and L5 ipsilateral dorsal roots could be visually identified by their size and, after opening the dura and applying a few drops of lidocaine to the cord surface, the ipsilateral dorsal roots from L1 to L6 were severed. We usually also cut a dorsal root rostral and caudal to the identified roots in order to ensure full deafferentation of the limb. The deafferentation was confirmed during the experiments by observing a lack of response to a strong pinch of the ipsilateral foot while observing a response to a similar pinch applied to the contralateral foot used to assess anesthetic state.

The animal was then transferred to the experimental setup and positioned (Fig. 1A). The posts cemented to the pelvis were attached to magnetic stands anchored to the table surface so that the pelvis was approximately level in medial–lateral and anterior–posterior directions. The threaded rod in the distal tibia was attached to a six-axis force transducer (Mini40; ATI, Apex, NC, USA). The transducer was mounted on a series of stages, which allowed it to be moved along the dorsoventral, mediolateral and rostrocaudal axes of the rat. Once the limb was attached to the transducer the only possible movement at the ankle was a rotation around the mediolateral axis. This rotation could be prevented or permitted by set screws. The ankle was positioned so that it was approximately in the same mediolateral plane as the hip center. The temperature of the rat was monitored with a rectal thermometer and maintained between 37 and 39°C using a heat lamp placed above the animal. Fluids were given to the animal immediately after the initial surgical preparation (dextrose 3 ml subcutaneously) and in between collections of force field data sets (approximately every 3 h).

We first determined recruitment curves for individual muscles, varying the current level and measuring the change in evoked force. Stimulation trains (75 Hz, 0.2 ms pulses, 500–1000 ms duration) were applied through the implanted cuff electrodes in order to evoke a muscle contraction with the limb fixed at a specific configuration in the workspace. The forces and moments measured by the transducer were collected (2000 Hz) using custom software in Labview for a period of 1 s before and 3 s after the stimulation train onset.

Stimulation trains were applied at a fixed interval of 1 min. Current was increased gradually from the threshold response until the evoked peak magnitude reached a plateau level. We measured the response for several current levels above the plateau level as well, making sure that the evoked force was not altered, in order to ensure lack of current spread to adjacent muscles. From these recruitment curves, we chose a current level that evoked a maximal activation of the muscle (between 0.06 and 0.45 mA). Once the current level was established, we then measured how the evoked forces and moments varied across the hindlimb workspace. The limb was fixed in a configuration, the muscle stimulated and the evoked forces and moments measured at the ankle with the attached force transducer. Two types of trials were collected for each muscle at a single limb configuration. In one trial, all degrees of freedom at the ankle were restrained (by securing the set screws; referred to as ‘restrained’ trials). In the other trial, rotation about the mediolateral axis of the ankle was permitted (by releasing the set screws; referred to as ‘released’ trials). By visual inspection, any such rotations were generally very small during evoked responses. Stimulation trains to any muscle were applied at an interval of 1 min. After collecting the evoked responses at a single limb configuration, the ankle was moved to a new position. We took care not to move the limb to the extremes of the workspace to avoid kinematic singularities or overstretching muscles. Stimulation was then applied to each muscle at the new position, for both restrained and released conditions. For each animal, 19–24 limb configurations were measured in order to span a range of limb configurations throughout the hindlimb workspace. Periodically through the collection of these responses the limb was returned to the original configuration and stimulation repeated to evaluate the effects of fatigue on the evoked responses. At the end of the experiment, the animal was killed by pentobarbital overdose and the position of the hip center of rotation was measured and the limb was dissected to find the distance from the hip to the knee center of rotation and the distance between the knee center of rotation and the rod attached to the foot. The mediolateral distance from the origin of the force transducer to the center of the foot was also measured for each animal in order to set the tool transform for the collected forces and moments.

We examined only the active force produced by each muscle, obtained by subtracting the averaged forces and moments measured prior to the onset of each stimulation train. We then found the peak force produced by each muscle and extracted the six-dimensional vector of peak forces and moments at each position. The set of these six-dimensional vectors measured at each limb configuration comprises a ‘force–moment field’ similar to the force fields that have been used previously (Giszter et al., 1993; Loeb et al., 2000). As described above, these fields are the result of the intrinsic force-generating properties of each muscle acting through the attachment of the muscle to the skeleton. In the approach described here we attempted to find a model of muscle action that is capable of explaining the observed muscle actions.

As shown in Fig. 1B, we modeled the rat hindlimb as a two-link serial chain consisting of the femur and tibia–fibula. The foot was omitted from analyses because we were concerned with muscles acting on the proximal joints. The hip was modeled as a ball joint and the knee was modeled as a universal joint with motion in flexion–extension in the sagittal plane and motion in pronation–supination of the tibia–fibula. There were thus five degrees of freedom in the limb. Constant centers of rotation were used for both the hip and knee joints.

If one considers the configuration of these experiments, it can be seen that there is a mismatch in the number of degrees of freedom within the hindlimb (five) and the measured forces and moments (six). Because of this, when the ankle was completely restrained by fixing the rotation around the mediolateral axis, the limb was overconstrained by the attachment to the force transducer. In such a situation, translating between intrinsic forces produced by the muscle and the forces and moments measured by the force transducer is not a well-posed problem, i.e. in general, there is not a unique mapping between the two sets of forces although the physical system has a unique mapping, dictated by the intrinsic stiffness of the force transducer and hindlimb (Mussa-Ivaldi and Hogan, 1991).

The M–P pseudoinverse is similar to the standard matrix inverse in that gives τ back, as can be easily seen from Eqn 5. Of the many possible values of ψ_a, Eqn 5 picks the one with the smallest magnitude: i.e. the least squares inverse of the matrix. However, this solution does not necessarily correspond to the true mechanical solution determined by the stiffnesses of the system and other pseudoinverses are also possible (Mussa-Ivaldi and Hogan, 1991). One important issue in considering M–P pseudoinverses is that the solution satisfying the least squares criterion depends on the choice of units used in the measurements. In our case, ψ_a consists of both forces and moments, which have different units. This is a well-known issue in robotics [e.g. appendix A, section 3.2 in Murray et al. (Murray et al., 1994)], for which there is a simple and standard solution, which is to choose an appropriate length scale to compare rotational quantities and translational quantities. In our case, the length scale is based on the size of the limb. We carefully chose units (newtons for forces and newton centimeters for moments) to make sure that both types of properties have similar orders of magnitude. Further, we verified that the results were minimally affected by changing the relative weight between two units within a reasonable range. However, because of these issues of the indeterminacy of the pseudoinverse, we compared this approach with a complementary one in which a pseudoinverse was not necessary.

Releasing this degree of freedom makes moments around this axis to be zero, but, in principle, it may cause the forces and moments measured along the remaining restrained degrees of freedom to differ from those measured with the limb fully restrained. However, we found that the forces and moments were only minimally altered when this degree of freedom was released (see Results).

A concern in fitting any model is whether the identified model overfits the data, explaining noise in the data rather than the structural variation. We therefore performed an n-fold cross validation. Each force field was divided into five sets of positions, chosen randomly. One set of positions was set aside and the remaining four sets were used to find model parameters, using the optimizations described above. We then evaluated the prediction of the identified model over the set of positions which were set aside. This procedure was repeated for each of the five sets of positions.

To quantify the quality of the fit to the observed data, we used a centered R² computed as one minus the ratio between the residual error variance in the model and the total variance of the data. An R² value near one indicates that there is minimal residual error, whereas an R² near zero indicates that the model performs the same as the mean of the data. Note that this is a conservative measure of fit quality in the present case, because predicting the mean force produced by the muscle is non-trivial but is a result of the specific origins and insertions found by the optimization.

Johnson et al. (Johnson et al., 2008) report a comprehensive musculoskeletal model for the rat hindlimb based on detailed anatomical measurements. We used these data plus the muscle architecture data of Eng et al. as the basis of a forward model to compute the force–moment fields in the present study (Eng et al., 2008). Johnson et al.'s variable joint center data, femur and tibia lengths, and muscle origin and insertion locations were incorporated in the model. This model will be referred to in this manuscript as the Johnson parameter model (Johnson et al., 2008). [Note that Johnson et al. have updated the model in SimTK, which is available online at the SimTK website (https://simtk.org/home/rat_hlimb_model) to include wrapping surfaces, made slight changes in some parameters, and the most recent model version considers the knee to move in only flexion extension. In the predictions shown here we used the model parameters and structure as described in Johnson et al. (Johnson et al., 2008).] As with our method above, the model has five degrees of freedom: X, Y, Z rotations for the hip joint and U, W rotations for the knee joint (see Fig. 1). To approximate the experimental conditions, we assumed that the knee was in the same mediolateral plane as the ankle (see the figures). Using the femur and tibia lengths from Johnson et al., we found the corresponding joint angles, allowing for the variable joint centers described by Johnson et al. These five angles were used to position the skeleton and determine the origin and insertion of the BFp and the gracilis posticus (GRp). Ankle force and moments were calculated using Eqn 4. J_a and J_m were determined numerically from the Johnson parameters for each set of joint angles. Because in the Johnson et al. manuscript joint centers are a function of joint angle, J_a and J_m as calculated here are slightly different from the analytic calculations used above. The magnitude of f_m was determined using a force–length relationship based on the data of Eng et al. (Eng et al., 2008). The shape of the function was approximated using a idealized sarcomere force–length function (Gordon et al., 1966) corrected for actin length in the rat (Rassier et al., 1999). The width of the function was scaled by muscle fiber length. The optimal length was set to muscle length. The maximal force was based on physiological cross-sectional area multiplied by 22 N cm^–2.

Fig. 2A shows an example of the forces and moments evoked from stimulation of the biceps femoris at a single limb configuration. As seen in the trial, the evoked response quickly reached a plateau level, which was sustained for the duration of the stimulation train then returned to the prestimulation baseline levels. Fig. 2B shows an example of the variation in the evoked force magnitude with the current level used to activate the nerve. In this example, the magnitude reached a plateau at 0.13 mA. Further increases of current past this level did not result in increases in the evoked force magnitude, suggesting that the muscle was maximally active. Fig. 2C shows the force and moment vectors in the X–Y plane taken at the time of peak response for the trial shown in Fig. 2A. We have separated the plots of forces and moments in order to make them easier to see, although in the data analyses and development of muscle models described below they were combined in a single six-dimensional vector. We also show the force and moment vectors evoked on repeated trials during the collection of a single force field in Fig. 2C. The vectors were collected at intervals of 8, 12, 12 and 14 min. The similarity of these vectors demonstrates the consistency of muscle stimulation in these experiments, although there is also fatigue over time because of the many trials of stimulation used here. Taken together, the results shown in Fig. 2 (see also Figs 3 and 4) demonstrate the reliability of the experimental preparation and collected data to characterize muscle actions.

Fig. 3 shows the force and moment fields evoked by stimulation of the GRp for each of the five animals, with the ankle movement fully restrained. The equivalent fields for the BFp are shown in Fig. 4. GRp produced a set of forces that drove the limb caudally, dorsally and medially. Fig. 3 shows that the force and moment fields measured for the GRp were very consistent across animals. Fig. 5A shows an example of two force and moment fields measured for one animal, either with the ankle fully restrained or with rotation around the mediolateral axis released. As can be seen in the figure, the moments around the mediolateral axis were minimal when the rotation was permitted while other degrees of freedom were only minimally affected.

Although the forces and moments measured for the GRp were generally consistent with what might be expected from the anatomy of this muscle (see below), the responses measured for the BFp were less expected. Although Fig. 4 shows that the responses evoked by stimulation of the BFp were similar across animals, the evoked forces were directed caudally in extension with only a minimal dorsal component. Such a caudally directed force suggests that the BFp acted predominantly as a hip extensor with minimal knee flexion. Fig. 5B shows an example of force and moment fields collected for BFp with mediolateral rotation either released or restrained, showing that the evoked forces and moments were minimally affected by releasing this degree of freedom.

The force and moment fields illustrated in Figs 3 and 4 reflect the intrinsic force-generating properties of each muscle and the attachment of the muscle to the skeleton. We, therefore, attempted to use the observed muscle actions to develop a parameterized muscle model for the GRp and BFp. Each model was described by eight parameters: the optimal muscle length, the maximal force, and the origin and insertion (each a point in three-dimensional space) of each muscle. The parameters that minimized the error between the observed and predicted forces–moment fields were found by optimization. We applied this optimization to the data obtained with the ankle either fully restrained or released, either using the Jacobian with M–P pseudoinverse (restrained analysis) or using the invertible Jacobian without z-moments (released analysis).

Fig. 6 shows the predictions of the model found for the GRp using the ‘restrained’ method for one animal. In each plot, the thick dark line indicates the predicted path of the muscle from its origin to its insertion. As can be seen in the plots, the identified model does a very good job of capturing both the forces and moments produced by stimulation of the GRp. The one exception to this good fit is for the z-moments: the observed response produced a negative z-moment, whereas the model predicted a positive z-moment (see Discussion). Because of this deviation in the z-moment, the R² value calculated for this case was quite low (R²=0.62). When this moment was excluded from the calculation, the R² value was much higher (R²=0.92). The origins and insertions identified by the model are also consistent with the expected anatomy of the GRp, with an origin near the midline of the pelvis and caudal to the hip, and with an insertion on the tibia partially towards the ankle. Fig. 8A (results from the restrained method are shown by empty symbols) shows the parameters for the origins and insertions found for the GRp in each animal. As can be seen in the figure, these parameters were very similar across animals (s.d.<0.15 mm), as might be expected based on the similarity of the force–moment fields shown in Fig. 3. Fig. 8C shows the predicted force–length functions for each animal, suggesting that in each animal the GRp was acting on the ascending portion of its force–length function.

Fig. 7 shows the results of the model found for the BFp using the ‘restrained’ method for the same animal shown in Fig. 3. As was the case for the GRp, the identified model did a very good job of capturing the observed actions of the muscle for both forces and moments. The origin of the model is also consistent with the expected anatomy of the muscle; it arises from the ischium caudal to the hip. Note the more lateral and dorsal origin of the BFp compared with the GRp. The insertion identified in the model was very close to the knee, consistent with the mainly hip extensor action suggested by the plots in Fig. 4. Fig. 8B (empty symbols) shows the origins and insertions found for the BFp model for each animal, again showing the consistency of the estimated parameters across animals (s.d.<0.17 mm). Fig. 8D shows that, similar to GRp, the model predicts that the BFp acts along the ascending portion of its force length function.

As described above, we also identified models based on the data sets with ‘released’ trials, in which the rotation around the mediolateral axis was permitted. These models were found using the invertible form of the Jacobian. Fig. 8A,B (filled symbols) shows the parameters found using these models from the released analysis for the GRp and BFp, respectively. Comparison of empty and filled symbols in plots shows that, in general, the model parameters found using the two approaches were very similar to one another.

The models using the ‘released’ method consistently had high R² values for both the BFp (0.88±0.02; mean ± s.d.) and the GRp (0.89±0.04), consistent with the good fit quality seen in the plots. The models found using the ‘restrained’ method for the BFp also had high R² values (0.84±0.06). For the GRp restrained analysis, the R² values were considerably lower (0.27±0.29), primarily because of the consistent offset in the z-moment as illustrated in Fig. 3. The R² values for the GRp increased considerably if this z-moment was ignored (0.85±0.06). Note that this is a centered R² measure, which considers only variations around the mean value in each force–moment dimension. Table 1 summarizes R² values for all five animals across the different conditions. Finally, we examined whether the models identified above were overfitting the observed data. We performed an n-fold cross validation procedure, fitting the data on 80% of the data then evaluating the fit on the remaining 20% of the data. We found that the quality of the fits was only minimally reduced when evaluating such cross validation fits, considering either the force–moment fields predicted using the ‘released’ method for the GRp (0.85±0.03) or BFp (0.85±0.06) or using the ‘restrained’ method for the GRp (0.81±0.06, ignoring the z-moments) and the BFp (0.85±0.03, ignoring the z-moments). This lack of overfitting suggests that the six-dimensional data vectors measured at multiple positions across the workspace place strong constraints on the eight parameters used in the model.

Taken together, the results show that simple muscle models consisting of point origins and point insertions are able to capture the actions of individual muscles in the rat hindlimb very well using either the ‘restrained’ or the ‘released’ procedures. Furthermore, the models identified by the two procedures and identified across animals were very similar to one another.

Fig. 9 shows the force and moment fields predicted from the Johnson parameter model (Johnson et al., 2008). As can be seen from the figures, the fields predicted for the GRp from the published parameters were very similar to those found in the present study by direct stimulation of the GRp (Fig. 3). However, the force–moment field predicted for the BFp is very different from the one observed experimentally using direct stimulation (Fig. 4). In particular, there is a much more prominent knee flexion action predicted for the BFp using the published parameters. This knee flexion action is a consequence of the insertion of the BFp being placed distally along the tibia, reflecting the anatomically distributed insertion of the BFp across the tibia. The results of this comparison suggest that the mechanical action of the BFp as measured here is better predicted by a more proximal insertion, such as those shown in Fig. 4 based on direct in situ measurements of BFp action, than by a more distal insertion, based on measurements of BFp anatomy.

As has been suggested in previous work (Giszter et al., 1993; Loeb et al., 2000), we found that the force and moment fields measured here give a succinct representation of muscle actions and their potential contributions to motor control. For instance, although standard characterizations might describe both the GRp and BFp as hip extensor and knee flexors (Greene, 1935), the force and moment fields shown in Figs 3 and 4 illustrate that the actions of these muscles are clearly distinct, not only in terms of their actions in the sagittal plane but also along other degrees of freedom. The description here involving six-dimensional forces and moments provides an especially rich description of muscle actions. It is likely that similar characterizations for other hindlimb muscles (Lawrence et al., 1993; Young et al., 1983) can help to highlight the unique contributions of individual muscles and provide insights into the coordination strategies implemented by the nervous system. The force and moment fields measured in response to stimulation of the GRp were consistent with the actions expected from anatomical inspection of this muscle and were similar to those predicted on the basis of parameters taken from a previous study using anatomical dissections (Johnson et al., 2008). This consistency helps validate the force and moment fields found here and the general approach taken in this study, in addition to helping validate the findings of that previous study. The fields measured in response to the BFp stimulation, in contrast, were more surprising. They reflected a predominant hip extension action, with a relatively small knee flexion contribution. This dominant hip extension action was also confirmed by the proximal insertion of the BFp estimated by our optimization procedures. Straightforward anatomical inspection of BFp suggests that it should produce a more significant knee flexion action than that measured here (Johnson et al., 2008). In fact, the fields predicted by parameters estimated in a previous study, using anatomical inspection, were clearly different from the ones measured here, reflecting the more distal insertion of the BFp estimated in that study. It is unlikely that the mainly hip extension action of the BFp observed here reflected an incomplete activation of all regions of the BFp because with increasing stimulation strength the evoked force magnitude reached a constant level, which did not change with further increases in strength (Fig. 2B). This difference between the insertion of the BFp estimated from anatomical dissection and the insertion estimated from measuring its action illustrates the utility of the approach described here in helping validate and refine musculoskeletal models.

Closer inspection of the structure of the BFp might provide an explanation for how this predominantly hip extensor action is consistent with the anatomy of BFp. Although the insertion of the BFp appears uniformly distributed along the tibia from the knee towards the ankle, reflecting the muscle to expose the medial surface reveals a well defined tendon near the knee into which a large portion of the muscle fibers in the BFp insert. The BFp is also thickest proximally and tapers more distally. Both of these features make it likely that this proximal insertion dominates the action of BFp when the entire muscle is stimulated. This division between a proximal region of the BFp with a well-defined tendon and a more distal region of the muscle with a distributed insertion along the tibia had apparently been recognized by Greene (Greene, 1935), who labeled the more proximal region the biceps femoris posterior and the more diffuse distal region the biceps femoris accessory. However, the exact division between these two regions is not obvious from simple inspection (Nicolopoulos-Stournaras and Iles, 1984; Parsons, 1896). In the cat, the BFp has been suggested to be divided into three compartments defined by separate innervation territories (Carrasco and English, 1999; Chanaud et al., 1991) and the divisions described here in the rat might similarly correspond to different compartments. Further experiments stimulating individual nerve branches to the BFp or partially cutting its insertion might help clarify the actions of different regions of the BFp.

It should also be noted that in order to isolate the muscles stimulated here we had to perform extensive dissections of the hindlimb. These dissections obviously interrupted connective tissue interactions between muscles, for instance between the BFp and the semitendinosus, the BFa and BFp, and the GRa and GRp, all of which have close adhesions. Similarly, in order to place bone screws on the distal tibia we interrupted the fascial sheath of the lateral limb. The importance of these intermuscular connections is not fully understood. In some circumstances, these interactions have been shown to significantly affect the forces produced by stimulation of individual muscles (Huijing and Baan, 2001; Hyde et al., 1999). During more physiological conditions, however, they have not been found to make a substantial difference in the force measured (Maas and Sandercock, 2008; Sandercock, 2005). In addition, some fascial sheets may need to be tensed by active muscles in order to transmit force normally. To characterize such effects, individual muscle actions might have to be studied on top of the background activity of other muscles. Future experiments might measure the force and moment fields with and without such connective tissue interactions or with and without other muscles pre-activated to gain a better understanding of the consequences of these nontraditional pathways of force transmission.

Beyond its use in validating and refining musculoskeletal models, we also demonstrated that the approach described here can be used directly to specify the parameters of a muscle model. This approach is related to previous work in which moment arms were estimated at single postures based on forces applied to tendons in cadavers (Lee et al., 2008; Pearlman et al., 2004; Santos and Valero-Cuevas, 2003; Valero-Cuevas et al., 2003; Valero-Cuevas et al., 2000) and in which planar force fields were used to find functional approximations of moment arm variations across the workspace (Loeb et al., 2000). The present study extends that work to use the forces and moments to specify a simple anatomical model of muscle action. We showed that a simple model of muscle action, based on a single origin and insertion along with the optimal length and force of the muscle, was able to describe the actions of both the GRp and the BFp very well. The good fits for the GRp with this model were consistent with the relatively circumscribed origin and insertion of this muscle, suggesting the validity of the approach used here.

However, we were surprised by the ability of this simple model, with a single origin and insertion, to explain the actions of the complex muscle BFp. We expected this muscle to show the limitations of this simple model structure, showing how the measurements from direct muscle stimulation could be used to justify the creation of distinct subregions and to specify their placement. As discussed above, it might be that this good fit reflected the dominant action of the more proximal tendon of the BFp, allowing it to be well described with a single insertion. Whether other complex muscles can be as well described with a single origin and insertion will require further experiments.

Ultimately, whether a complex muscle such as the BFp should be represented by a single insertion or multiple insertions and where those insertions should be placed depends not only on being able to explain the mechanical action of the muscle but also on being able to explain the activation of the muscle by the nervous system. For example, even though the present experiments suggest that we could explain the mechanical action of the entirety of the BFp by a single insertion point, it is probable that the nervous system does not activate the BFp homogenously. Previous work in the adult cat (Chanaud and Macpherson, 1991; Pratt et al., 1991) and our own work in the neonatal rat (Klein et al., 2010) have shown that regions of the BFp are activated differentially in behaviors such as locomotion and posture. Obviously, the mechanical consequences of this differential activation cannot be explained using a musculoskeletal model of BFp with a single origin and insertion. Such considerations emphasize the iterative process of refinement and elaboration involved in creating musculoskeletal models, requiring a combination of both biomechanical and neurophysiological studies.

The results of this study were mainly focused on estimating the origins and insertions of the GRp and the BFp, based on the structure of the measured force–moment fields. These origins and insertions determine how the direction of the six-dimensional vector measured at the ankle varies across limb configurations, which is a crucial aspect of all musculoskeletal models. Our models also included estimates of the force–length properties of each muscle, as defined by the maximal force and optimal muscle length parameters. These properties influence the magnitude of the vector at each limb configuration but do not affect the direction of the predicted vector. The model of force generation used in this study only included parameters defining the maximal force and optimal length of each muscle and did not include a parameter for the resting tendon length (Zajac, 1989). Because we only made isometric measurements in these experiments and because these muscles are likely to be acting mainly on the ascending portion of their force–length functions (see Fig. 8, and from predictions from Johnson model parameters), including this extra parameter is not expected to increase the fit quality of the identified models, consistent with the very good fits shown here. For other muscles, such as distal muscles with long tendons or muscles operating on both the ascending and descending portion of the force–length function, it should be possible to estimate the tendon length and increase the quality of the observed fits by including this parameter. Note, however, that because of the independence between force–moment direction and intrinsic force generation described above, this simplification would not be expected to alter the estimated origin and insertion of each muscle. Similarly, although equations other than the one used here could be used to characterize the force–length relationship of muscle, it is unlikely that their use would affect the estimated anatomical parameters. Further experiments in which muscle actions are characterized dynamically or across a wider range of conditions might make it possible to estimate these force-generating properties of muscles more completely.

Note also that although the experiments described here involved extensive dissections to isolate individual muscles, other less invasive methods of estimating muscle actions have been developed (Kutch et al., 2010). Such less invasive methods might be used in combination with the computational analyses described in this study to construct musculoskeletal models.

The present study demonstrates the potential utility of in situ measurements of muscle actions to validate and develop musculoskeletal muscle models. However, it is important to point out several limitations of the present study. First, the measurements made here are obviously subject to experimental error. This error arises for a number of reasons, including imperfect mechanical isolation of the hip, muscle fatigue both within and between trials, and idealized kinematic models of joints. For instance, the z-moments observed for the GRp could not be explained with the models considered here. In fact, the negative z-moments observed from stimulation should be impossible for any muscle if the hip was completely mechanically fixed. Obtaining ideal mechanical fixation is clearly difficult because of the considerable internal flexibility in the sacrum and pelvis: despite placing bone screws on the contralateral pelvis and on the rostral ipsilateral pelvis, there could still be some movement of the hip during muscle stimulation because of internal flexion of the pelvis. We were unable to restrain the skeleton further without interrupting the insertion of the BFp on the pelvis. Although this imperfect fixation can apparently alter the measured actions of muscles, it is important to note that these effects did not substantially affect the models identified in these experiments: if we excluded the z-moment from the optimizations, the anatomical parameters estimated were affected only minimally. Using more accurate joint models and tracking the exact skeletal configuration at each limb configuration would also potentially increase the accuracy of models developed using the methods described here.

The ‘top down’ approach described here to estimate the static properties of muscles might, in principle, also be used to estimate the dynamic properties of muscles, taking measurements of the behaviour of the intact muscle when the limb is being moved through a variety of trajectories, then using those measurements to estimate dynamic properties of muscles. Whether the increased complexity of such experiments would still provide useful estimates of these properties or whether straightforward morphometric measurements of muscles would be simpler and more accurate is not clear at present. Regardless of whether such experiments might be useful for estimating musculoskeletal models, it is clear that they would be useful for validating models created from other procedures. Evaluating the utility of this approach for dynamic properties of muscles will require further experiments.

It will be interesting in future work to evaluate the sensitivity to these sources of error of the models identified using the top down approach of this study. Whether or not the consequences of these measurement errors are greater in this top down approach than they are in more traditional ‘bottom-up’ approaches based on anatomical measurements (Santos et al., 2009), it seems clear that combining the information from both approaches can be a useful way to create and refine complex musculoskeletal models.

 
• BFp

  posterior head of biceps femoris muscle

 
• d

  the eight-dimensional vector of muscle model parameters

 
• F_m

  (3×1) vector of forces

 
• f_max

  maximal force

 
• GRp

  gracilis posticus muscle

 
• J_a

  (6×5) Jacobian relating muscle force–moments measured at the ankle to joint torques

 
• J_a

  (5×5) Jacobian for the ‘released’ method

 
• J_m

  (3×5) Jacobian relating muscle force to joint torques

 
• l

  muscle length

 
• l_o

  optimal muscle length

 
• τa

  (5×1) vector of joint torques needed to produce ψ_a at the endpoint

 
• τm

  (5×1) vector of joint torques created by the musce force

 
• ψa

  (6×1) vector of forces and moments measured at the ankle

 
• ψa

  (5×1) reduced vector of forces and moments for the ‘released’ method

 
• ψ_i(d)

  ankle force and moments predicted for the ith configuration