Method for determining musculotendon parameters in subject-specific musculoskeletal models of children developed from MRI data

Accurate knowledge of muscle-tendon parameters in biomechanical models is critical for accurate simulation and analyses of human movement. An excellent example of this is the creation of subject-specific models from magnetic resonance imaging (MRI). When Hill-type muscle models are used to calculate muscle forces, the determination of muscle attachment points, optimal fiber length, tendon slack length and maximum isometric force all have a significant influence on the joint moment-angle behavior of the model. In the present study a method was developed for customizing the values of muscle-tendon parameters in a generic model to create subject-specific biomechanical models from MRI. The method was applied by generating musculoskeletal models for the biomechanical simulation platform OpenSim, but the workflow is equally well applicable to other simulation platforms. New computational algorithms are described for identifying joint centers and for reconstructing the centroids of the muscle bellies from MRI. A process is also described for the extraction of the muscle paths and for identifying the positions of ‘via-points’ used to model muscles wrapping over bones. Finally, a new algorithm is described for adjusting the values of optimal fiber length, tendon slack length and maximum isometric force based on a comparison of the model results with experiment. We tested our computational algorithms by developing subject-specific biomechanical models of five typically developed children (age 9.5±1.7 years) from MRI. The joint moment-angle relationships calculated for the subject-specific models were similar to those determined for corresponding scaled generic models. The results indicate that the proposed methodology is suitable for developing subject-specific models of healthy children. Future studies should investigate how abnormalities of the musculoskeletal system, such as tibial torsion and muscle spasticity, can be integrated into the modeling process.


Introduction
Musculoskeletal modeling and simulation is used to study human movement and to better understand the biomechanical principles of human gait [1]. Generic models of the musculoskeletal system have been developed from anatomical and biomechanical studies of cadaveric specimens [2,3]. Although these models have been applied in many previous studies, skeletal and musculotendon parameters are different for each individual, resulting in the unique capability of each muscle to generate joint moment profiles. Gait-marker-scaled generic musculoskeletal models can be used to estimate a subject's musculoskeletal geometry, but this approach does not account for subject-specific geometry [4,5]. Subject-specific models need to be created for the assessment of individual musculoskeletal geometry and analysis of muscle function during motion.
Modeling based on magnetic resonance imaging (MRI) data can integrate subjectspecific details into musculoskeletal models [3,6,7]. Bone structures can be identified from MR images for more accurate modeling of body segments, their connecting joints and wrapping surfaces that define the path of a muscle spanning a joint. Muscle and tendon structures as well as their attachment sites can be identified and muscle moment arms derived. In order to accurately simulate the joint moments generated by a muscle, geometric parameters of the muscle-tendon actuators that influence its force generation must be determined. In the commonly used Hill-type muscle model [8], these parameters are tendon rest length, optimal muscle-fiber length, physiological cross section area (PCSA) and pennation angle [9,10]. Determination of tendon rest length and optimal muscle-fiber length from MRI images is challenging due to the difficulty of distinguishing the aponeurotic part of tendon from the muscle belly, and because the muscle is usually not at its optimal length during imaging. Because the resulting muscle behavior is quite sensitive to parameter changes [11], it is important to select these parameters carefully. It has been shown that tendon rest length and optimal muscle-fiber length do not scale linearly with bone length [12].
Maximum isometric muscle force is a function of muscle length and hence joint angle. Winby et al. [13] evaluated different approaches to subject-specific scaling of musculotendon parameters and suggested preserving the muscle's normalized fiber length (fiber length divided by optimal fiber length) between the unscaled and scaled models over the joint range of motion.
Scheys et al. [14] published an atlas-based method for defining muscle paths from MRI. One limitation of their approach is that it does not take into account the possibility that the subject's joints are not uniformly in the neutral position (e.g., the hip and knee fully extended and the ankle in the neutral standing position), which is often the case for subjects who have musculoskeletal disorders such as children with cerebral palsy.
In this study, we introduce a method to scale a generic model of the skeleton and joints, in addition to the lines-of-action of the muscles and their associated parameters, so that subjectspecific musculoskeletal models can be built in a systematic and reproducible manner. We used the method to generate subject-specific models of typically developed children based on a generic model of an adult. The generic model was scaled with subject-specific adjusted parameters derived from MRI. A semi-automated workflow was developed to extract the muscle paths from manually segmented MR images of the test subject and to specify the values of all required Hill-type muscle-tendon parameters appropriately.

Methods
The following sections describe the process of customizing a generic biomechanical model of an adult to obtain subject-specific models of children based on MRI data. We implemented the method in Matlab (1994 Mathworks, Inc., Natick, MA, USA) and used the model structure of the biomechanical simulation software platform OpenSim [15].

Selection and preparation of a generic biomechanical model
The generic biomechanical model of an adult was based on a model published by Arnold [16], which was developed on the basis of data obtained from 21 cadaver specimens. The model includes geometric representations of the bones, kinematic descriptions of the joints, and Hill-type models of 88 muscle-tendon compartments. We reduced the amount of actuators to a set of 46 muscles by lumping some muscles into single actuators. This was done because the individual compartments of some muscles were difficult to identify separately in the MRI dataset when segmenting the data manually. For example, the three compartments of gluteus maximus were combined into a single muscle. The muscle attachment sites of the new combined muscle as well as the optimal muscle-fiber length and tendon slack length were manually adjusted to obtain a similar joint moment-angle behavior displayed by the original muscle group.
For each muscle, muscle-fiber length (l g m.z ) and tendon slack length (l g t.z ) were extracted with the generic model (superscript g ) positioned in an erect standing position (i.e., with the hip, knee and ankle joints positioned at 0 degrees; subscript .z for zero degrees).

Subject data
For the subject-specific models we collected MRI data of five typically developed children (age 9.5 ± 1.7 years, height 1.34 ± 0.075 m and weight 30.3 ± 3.8 kg). The study was approved by the Local Human Ethics Board of Vienna and the parents of the children gave informed consent for participation.

MR-imaging and segmentation
MR images were collected at the MR/CT Institute Schmidt GmbH&Co KEG using a Siemens Symphony Maestro Class 1.5 Tesla scanner. The scanning parameters were selected to enhance the brightness of fatty tissue in order to make the boundaries around muscles more visible. Voxel sizes of 0.9 × 0.9 × 5.0(3.0) mm for the axial images and 1.8 × 1.8 × 8.0 mm for the sagittal images were used, which are comparable to the voxel sizes used by Spoor and van Leeuwen [17] and Arnold et al. [6]. The subject was placed in the prone position in the MR scanner with the hips and knees fully extended and the ankles in the neutral (standing) position. Data were recorded from both legs.
To define the bone and muscle surface geometry, four series of transverse images were obtained using a body coil. In the first series, 3 mm slices were obtained from the lower lumbar vertebrae down to the lesser trochanter of the femur. Three series were then recorded using 5 mm slices along the shaft of the femur and tibia from the lesser trochanter down to the  calcaneus. Two series of sagittal images with 8 mm slice thickness were also obtained assist with identifying the muscle structures and attachment sites during manual segmentation. The commercial software package Amira (Visage Imaging) was used for the image preparation. Nitro capsules mounted on the skin during MRI served as markers to align the corresponding image series for stacking. The stacked images were then merged using an included function of Amira to obtain a single full dataset defining the muscle and bone surfaces of the pelvis and legs. The boundaries of the pelvis, femur, tibia, and of all the major muscle groups in the lower limb were identified manually in each individual MRI slice. This procedure was repeated for both legs. Voxels within each boundary were added to the specific muscle's volume point cloud and certain anatomical landmarks identified and marked.
The segmentation process resulted in a three-dimensional matrix with specific values for identifying each marked structure. The three values indicating the position of each voxel inside the matrix multiplied by the voxel size corresponded to the metric position of the voxel inside the point cloud volume. The data were imported into Matlab, where the dataset was automatically analyzed for preparing the biomechanical model.

Analysis of MRI data-extracting the anthropometry and scaling
Virtual markers were placed at a number of anatomical landmarks as indicated in Table 1. The coordinates of these markers were extracted and used to derive the joint centers of the hip, knee, and ankle joints as well as the center of the pelvis in global coordinates. Based on the global coordinates of the joint centers the orientation of each segment was calculated in relation to the adjacent proximal segment.
The virtual MRI marker set described in Table 1 served as the basis for defining the anthropometry of each subject's lower limbs. The ratio of the MRI marker distances to the corresponding marker distances defined in the generic model was used to scale the adult generic model to the subject-specific model of each child. In a previous study [5] we showed that the scaling result based on a virtual MRI marker set is more accurate than scaling via skin mounted gait markers.

Analysis of MRI data-correction of limb orientations in MRI
We defined the center of the pelvis as the origin of the global coordinate system for the MRI data. The medial/lateral axis was oriented to the right ASIS and the anterior/posterior axis to the midpoint of the two PSIS markers. All coordinates defining the structure of the model skeleton were defined in the local coordinate system of the body segment in question. The origin of each body segment was set to the proximal joint center where the segment is connected to the next proximal segment. For example, all coordinates on the thigh (reference body femur) are given in a coordinate system placed and pivoted in the hip joint center.
Subjects lying in an MRI machine do not always have their joints exactly in the neutral position (i.e., 0 degrees for all joint angles). To extract the locations of the joints, muscle attachment points and muscle paths from the segmented MRI dataset, we took the limb positions into account. We calculated transformation matrices describing orientations and positions of the reference coordinate systems with respect to each other. These transformation matrices were referenced to a kinematic chain, and so it was possible to virtually retransform all joint angles to a neutral position. In this way, all desired locations from the MRI were obtained both in a global reference coordinate system and the local coordinate system of the segment in question.

Analysis of MRI data-extracting muscle geometry
To obtain a muscle's path, the central line throughout the segmented volume of each muscle was calculated by connecting the geometric centers of the slices and smoothing the resulting line with a running average filter using a window size of three slices.
As described in Sect. 2.5 above, the joint center locations and the muscle centerlines were determined for the neutral standing position. We reproduced the muscle paths by defining all points proximal to a joint in the coordinate system of the proximal body, whereas all points distal to a joint were first expressed in the slightly twisted coordinate system of the distal body segment and then transformed to the neutral standing position.
The volume of each muscle (V m ) was calculated by multiplying the summed number of voxels for each muscle by the voxel size used during MRI recording.

Muscle path modeling
The origin and insertion sites of each muscle were defined in the model using the centroid of the muscle derived from the MR images. For each muscle the most proximal point of the centroid was defined as the origin and the most distal point as the insertion. The subjectspecific muscle attachment points were manually inspected to correct errors that occurred during segmentation of the MRI images.
Via-points [1,2] are commonly used in biomechanical simulations to model the manner in which a muscle wraps around other muscles and/or bones within a defined range of motion of the joint. If the joint angle occurred in the specified range, then the muscle's path was defined by an additional via-point (P via ). If, on the other hand, the joint angle occurred outside the specified joint range of motion, then a proximal point (P pro ) and a distal point (P dis ) of the modeled muscle's path were connected by a straight line. Because the coordinates of the proximal point, distal point and via-point can be given in different reference frames, the via-point coordinates were found as follows: First, the locations of the proximal and distal points of the muscle's path were transformed into a local coordinate system with its origin in the processed pivot joint center (J ) and with one axis matching the rotation axis of the joint, thus: The distal point (P J dis ) was then rotated about the pivot joint center (J ). We used the joint angle range specified in the generic model, as this resulted in smooth joint moment-angle The via-point was placed on L at the shortest distance to the joint center, J . To compute the shortest distance, let P via be the base of the perpendicular extending from J to L as shown in Fig. 1. The vector P pro to P via is the projection of the vector from P pro to J onto L, as shown in the diagram. The global coordinate of the via-point is then given by Alternatively, computational algorithms used in other simulation platforms can be applied to calculate the necessary muscle wrapping paths [18].

Calculation of geometric muscle parameters and maximum isometric force
When a model is scaled in OpenSim, muscle-fiber lengths and tendon slack lengths are scaled so that they remain the same percentage length of the scaled total actuator length. If a muscle's attachment sites and path are changed, the length of a muscle-tendon unit is altered, but the values of optimal muscle-fiber length and tendon slack length are not automatically altered. Changing the muscle's attachment sites will change the muscle's joint momentangle behavior because (1) the muscle's moment arms are altered, and (2) the muscle no longer operates on the same portion of its force-length curve. This results in an unwanted shift in the optimum joint angle for the muscle, i.e., the joint angle at which the muscle has its optimum length. In the approach taken below, the optimal muscle-fiber lengths and tendon slack lengths are chosen so that the optimal joint angle is identical with that given by the generic model. The magnitude and shape of the joint moment over the joint range of motion depends on the maximum isometric force of the muscle, which in turn is a function of the muscle volume (determined from MRI) and the moment arm of the muscle which is determined by the muscle path.
It was assumed that the optimal joint angle at which the muscle fibers have their optimal length is the same for an adult and a nine-year-old child. This assumption is based on the fact that the joint angle trajectories during gait are similar for adults and children who are nine years and older. Assuming the shape of the muscle force-joint angle curve is the same in a healthy child as an average adult, the relation of muscle-fiber length to optimal fiber length is preserved for both groups at any joint angle. This approach is similar to the method used by Winby et al. [13] who obtained the best results when they preserved the normalized muscle-fiber length over the whole joint range of motion. These researchers also showed that the results differed only slightly when two different postures were selected for evaluation. We therefore evaluated the muscle parameters with the subject-specific model positioned in the normal standing position and in the position that defined the optimal joint angle for each muscle.
The following steps were used to determine the muscle parameters for the subjectspecific model: -Obtain the optimal muscle-fiber length for the subject-specific (l s m.0 ) model directly from OpenSim, where all geometric muscle parameters are scaled in proportion to model anthropometry.
-Compute the total muscle tendon length (l s mt.z ) with all joint angles in the subject-specific model set to zero degrees (i.e., the reference position) -Calculate an elongation factor for each tendon (ε t ), defined as the ratio between tendon slack length and the tendon length calculated with the generic model placed in the reference position, thus: -Calculate an elongation factor for each muscle-fiber length (ε m ), defined as the ratio between the optimal muscle-fiber length and the muscle-fiber length calculated with the generic model placed in the reference position, thus: -Compute the muscle-fiber length for the subject-specific model, l s m.z , by multiplying the optimal muscle-fiber length in the subject-specific model by the elongation factor for each muscle, thus: -Compute the tendon slack length for the subject-specific model, l s t.0 , by multiplying the tendon slack length in the subject-specific model, l s t.z , by the elongation factor for each tendon. Note that the muscle pennation angle is also taken into account here, thus: Thus, which can be re-written as Maximum isometric muscle force can be obtained by dividing the muscle's volume, V m , by the optimal muscle-fiber length, l s m.0 . The value of maximum muscle stress, σ , was taken to be 33 N/cm 2 [10]. Thus, the maximum isometric force of each muscle is given by

Model validation
To validate the modeling method we compared the magnitudes and overall shapes of the joint moment-angle curves of the subject-specific models to results obtained from a scaled generic model as well as experimental joint moment-angle data (Figs. 4 and 5). The generic model was scaled to the average size of the subjects (∼75%) and the maximum isometric forces were adjusted to muscle strengths corresponding to a reduced body mass (56%). We calculated the maximum joint moments exerted by each muscle over the range of joint motion, as well as the passive joint moments of the antagonistic muscles. To calculate the net maximum joint moment, we added only the 'positive' contributing sections of the joint moment-angle curve of each active muscle and deducted the passive moment of the antagonistic muscles.
To validate our results, we used averaged data obtained from the isometric joint moment measurements of the five typically developed children for hip flexion, hip abduction/adduction, knee flexion/extension and ankle plantar/dorsiflexion. The results of both legs of each child were included and averaged in the analysis. For hip extension, we used data measured by Anderson [19] that was scaled to the magnitudes of the joint moments measured for the children using Eek's [20] method. Overall magnitudes can be compared to the data reported by Eek [20] for specific joint angles. The results for hip flexion/extension were measured and calculated in a supine position with the knee flexed 90 • , for abduction/adduction supine with straight leg, knee flexion/extension sitting with 60 • hip flexion and for plantar/dorsiflexion while sitting with hip and knee 60 • flexed.

Illustration of modeling steps
To illustrate the modeling steps in Figs. 3 and 2, the graphs labeled 'scaled generic muscles 1' show the resulting joint moments over the range of motion once the generic model Fig. 3 Hip-flexion moment-angle curve calculated in the model was scaled to the height of an average nine-year-old child (∼75%) and the maximum isometric forces were set to approximately ∼56% of the generic adults strength. In the graph titled 'scaled generic muscles 2', also the via-point was corrected according to the new muscle path after scaling and the tendon slack length was set to the correct value as described in Sect. 2.8 (see (9)).
The graph titled 'Specific Model 1' was obtained once the subject-specific muscle attachment points were adjusted to the data derived from the MR images and the maximum isometric muscle forces were then set to the subject-specific values obtained from the MRI data. The discontinuities observed at a hip flexion angle of 50 • in Fig. 3 and at a knee flexion angle of −15 • in Fig. 2 were caused by incorrectly positioned via-points. After adjusting the position of these via-points and correcting the geometric muscle parameters, the final joint moment-angle curve (Fig. 3, Specific model 2) can be seen to have a similar shape to the joint moment-angle curve computed in the generic model.
The measured subject-specific attachment sites for vasti (Fig. 3) lie about 10 cm more proximal than in the generic model after scaling. This results in an unusual joint moment profile for knee extension as shown in the graph titled 'Specific Model 1'. After adjusting the via-points and tendon slack lengths in 'Specific Model 2' a reasonable match was obtained to the shape of the joint moment-angle curve computed in the scaled generic model. The discontinuity observed at a knee flexion angle of −80 • is model-specific and results from an incorrectly positioned attachment site for the vasti on the patella, which was not corrected after scaling.

Comparison of joint moments over the range of motion
We calculated the joint moment-angle curves for the subject-specific models of five children (age 9.5 ± 1.7 years, height 1.34 ± 0.075 m and weight 30.3 ± 3.8 kg) using the 'Analyze' function available in OpenSim. In Figs. 4 and 5, the averaged results of the subject-specific models, which were developed from MRI data according to the proposed method, are compared to the results of a scaled generic model.  [19]. Dot with error bar: Result from Eek [20], manual dynamometry on children The results show that the magnitudes and overall shape of the joint moment-angle curves in the subject-specific and scaled generic models are very similar for the joints examined in this study. However, some differences can be observed between the simulated joint moments and experimental measurements.
Compared to the result of Eek [20], a reduction in hip flexion moment can be observed both in the model and our experimental data. Eek measured the data with a hand dynamometer with the other degrees of freedom of the measured leg not completely locked. To examine the effect of this limitation, we calculated the maximum hip flexion moment in the model with a slight hip adduction of 10 • that might have helped the subject to stabilize and a possible knee angle of greater than 90 • that shortens the hamstrings tendon and therefore reduces the hamstrings passive joint moment during hip flexion. With this configuration, the maxi-  [20], manual dynamometry on children mum joint moment computed in the model was up to 30% higher, particularly when the hip flexion angle was increased.
The generic model shows a decreasing ability of hip abduction the more the hip is adducted. This behavior is also reproduced by our modeling method. Contrary to this our experiments show increased hip abduction ability when the hip is adducted. The underlying generic model needs to be examined more closely and revised regarding this issue, especially in terms of muscle paths.
The measurement for knee extension was performed in a sitting position of the test person where the hip was flexed to approximately 60 • . In this posture, the knee flexors that also cross the hip joint (e.g. hamstrings) are lengthened, which leads to a rising resistive moment when extending the knee. Decreasing tendon slack length for the hamstrings by around 5 cm generated strong passive joint moments at smaller knee angles. This in turn reduced the net knee extension moment in a way that the shape of joint moment calculated in the model more closely resembled the result measured for the subjects. This issue will also need to be examined more closely and revised in the underlying generic model.
Differences in the magnitudes of the joint moments may reflect differences in the maximum isometric muscle forces assumed for the individual subjects as well as variations in muscle moment arms that result from geometric differences between the subject-specific models and the scaled generic adult model. For example, the large difference in the magnitude of the calculated dorsiflexion moment was based on the reduced maximum isometric forces assumed for the extensor digitorum and tibialis anterior in the MR-based models, which were about two-thirds of the values assumed in the scaled generic model.

Conclusions
Our approach to subject-specific biomechanical modeling has led to reasonable results, at least as far as reproducing the moment-angle behavior of joints is concerned. As children at the age of about eight years have similar gait patterns to adults, it was assumed that the lower-limb muscles of children and adults have their optimal lengths at the same joint angles. Maximum isometric muscle forces, tendon slack lengths and moment arms were based on the data derived from MR images. The resultant joint moments calculated in the subjectspecific models were similar to the moments determined from the scaled generic model, which seems reasonable because the gait patterns of the two subject cohorts were also similar. However, the functional roles of the individual muscles may be different, but this can only be established by performing further analyses of the model simulation results (e.g., by calculating induced accelerations of the joints and the body's center of mass). The results presented in this paper show that the proposed methodology is suitable for developing subject-specific models of healthy children. Future studies should investigate how abnormalities of the musculoskeletal system such as tibial torsion and muscle spasticity can be integrated into the modeling process.
While we have not applied the proposed methodology to other modeling platforms (aside from OpenSim) or to models of different anatomical structure, the approach we have taken can easily be applied to these scenarios as well. The computational algorithms are implemented in Matlab, and the implementation is based on the structure of OpenSim models, with model parameters changed via Matlab. If an interface to another software platform is available, then the algorithm could be easily applied to read and write model parameters using Matlab.
One limitation of this study relates to the manner in which muscle-tendon attachment sites were selected. It is challenging to confidently identify tendon attachments to the bones, since the attachment site is usually not a single point, but an extended area (e.g., the medial portion of gluteus medius and the medial portion of iliacus) or a finite length along the bone (e.g., vastus medialis and the medial portion of adductor magnus). The centroid of the muscle belly and tendon can be used to model a muscle's path. However, it remains unclear whether this muscle path is identical with the muscle's line of action. In particular, when the attachment site of the tendon on the bone is extended, it is difficult to determine which parts transfer muscle stress to the skeleton during movement. By selecting the most proximal and distal occurrences of a muscle-tendon structure as attachment sites and calculating the centroid as the line of action of the muscle, we obtained a rather simple but repeatable method for reproducing the line of action of a muscle.
After accurate segmentation, the customization of a model took only about five minutes using the proposed method. The time taken to generate subject-specific musculoskeletal models can be reduced even further by automating the steps needed to reproduce the path of each muscle. Nonetheless, most of the time taken to build a subject-specific musculoskeletal model is related to manual segmentation of the MR images. Atlas-based methods may help in this regard, but if subjects differ from the norm due to individual pathologies, then this approach may not assist in reducing the overall time required for model development.
Our future work involves applying the proposed method to build subject-specific models based on MRI data obtained from children with cerebral palsy and to use these models to perform simulations and analyses of these subjects' gait patterns.