A Cadaveric Comparison of the Kinematic and Anatomical Axes and Arthrokinematics of the Metatarsosesamoidal and First Metatarsophalangeal Joints

Essential to our daily ambulation and general foot function, the first metatarsal and proximal phalanges, along with the medial and lateral sesamoids, are responsible for transmitting a significant percentage of bodyweight during gait [1]. The sesamoids embedded in the flexor hallucis brevis tendon glide on the distal plantar metatarsal head in the metatarsosesamoidal joints (MTSJ). They simultaneously act to increase the moment arms of the flexor hallucis longus and brevis tendons of the first metatarsophalangeal joint (MTPJ1) and protect the metatarsal cartilage from compressive and shear stresses [2]. Dysfunction and degradation in these joints are of interest to clinicians studying and treating maladies such as hallux valgus, hallux rigidus, and clawed hallux [1,3–6]. By quantifying the motion of these small bones in both healthy and pathological joints, clinical treatment strategies and our understanding of bone function in the foot will evolve.

Precise data describing the motion of the sesamoids and first proximal phalanges are sparse, due to limitations in traditional, noninvasive techniques for assessing foot and ankle biomechanics. Their small sizes and locations encased in supporting tissue structures, make the sesamoids particularly challenging targets for conventional imaging methods; they require static, partial-weightbearing, tangential radiographic views [1,7] and they cannot be tracked using optical motion capture. As a result, MTPJ1 and MTSJ studies have been limited to cadaveric dissections that required disruption of the native anatomy [8–10] or low-load in vivo magnetic resonance studies of the plantar plate [11–13] to quantify bone motion. Various motion capture models have been developed to improve forefoot motion quantification [14–16], but have inherent limitations in the degrees-of-freedom that can be tracked accurately for each bone. Recently, weightbearing computed tomography (WBCT) technology has advanced with cone-beam sources and high-efficiency detectors, offering low-dose, high-resolution volumetric scans of clinical patients' extremities [17,18]. Such systems allow for assessing the full three-dimensional shapes of the foot bones and their six-degree-of-freedom kinematics via in vivo or controlled cadaveric experiments.

Given the current lack of precise forefoot kinematic data, there are no standardized methods for establishing coordinate systems and reporting kinematics for the MTSJ and MTPJ1 as in other joints, such as the knee and ankle [19–21]. Coordinate systems may utilize anatomical landmarks, the fitting of geometric primitives, or calculating the inertial axes (principal components) of the bone [22–25]. To mitigate crosstalk and between-subjects variances in the resultant kinematics that confound or dilute differences between cohorts, the natural variation in bone shapes and the role of articulations should be considered when defining anatomical coordinate systems. Such articular surface-based coordinate systems have been utilized in the upper extremity, knee, and hindfoot [26–30]. Additionally, motion-based finite helical axes, which have characterized knee and wrist joints [31,32], invite comparison to these coordinate systems based on osseous morphology. Quantitative comparisons in resultant kinematics between coordinate systems and their sensitivity to crosstalk have been performed in the knee joint [33–35]. A quantitative method is needed for assessing and justifying the definitions of anatomically based joint coordinate systems, especially as technology for foot joint tracking improves and facilitates more precise analyses. The primary aims of this study were to quantify the arthrokinematics (estimations of joint contact foci), instantaneous centers of rotation, and joint kinematics of the MTPJ1 and MTSJ, and to compare the orientation and distribution of the joint mean finite helical axes (MFHA) to the medial-lateral axes of first metatarsal coordinate systems based on two methods: principal components analysis (PCA) and fitting cylindrical primitives (CYL) to the bony anatomy. For MTPJ1 and MTSJ, their primary motions are flexion–extension about the distal metatarsal head, so the MFHA representing these motions should closely align with the coordinate system medial-lateral axes. PCA equally weights the entire metatarsal in the coordinate system definition, including proximal portions away from the distal joints of interest. A localized fitting method like CYL can adapt to variations in torsion of the distal metatarsal head relative to the shaft [24,36–38]. We hypothesize, first, (H1) that vector directions of the MFHA of MTPJ1 and MTSJ will exhibit stronger spatial correlations with medial-lateral metatarsal axis definitions that are based on the CYL method compared to PCA, and, second, (H2) that coordinate systems defined using the CYL method will produce MTPJ1 rotation data with less kinematic crosstalk and less intersubject variability compared to data derived from PCA-based coordinate systems.

Twelve fresh-frozen tibia-foot specimens were obtained from accredited tissue banks (Innoved Institute, LLC, Elk Grove Village, IL, and Lonetree Medical Donation, Centennial, CO) and screened for abnormal bony alignment or deformity via donor history, manual palpation, and radiographic imaging. Normal range of motion and lack of osteoarthritis was confirmed by an orthopedic surgeon using palpation, radiographs, and WBCT (pedCAT, CurveBeam, Hatfield, PA; isotropic voxel size of 0.3 mm). Two specimens having bipartite sesamoids were omitted. For the remaining ten specimens (age: 46.3±17.2 years, body mass index: 25.8±4.7, sex: six male, four female), all subsequent scans were performed with the pedCAT. Specimen-specific segmentations in mimics software (v21, Materialize, Leuven, Belgium) by a single operator produced a volumetric template of each specimen's first metatarsal, proximal phalange, medial sesamoid, and lateral sesamoid for tracking in subsequent scans.

To express kinematic and arthrokinematic data between specimens in a common frame of reference, standardized specimen-specific anatomical coordinate systems were defined for each bone. In the metatarsals, two methods of embedded coordinate systems were compared: one from the principal components of the entire bone surface [22,23] (PCA) and the other by fitting cylindrical primitives to the metatarsal surfaces [24,25] (CYL). For our CYL protocol (Fig. 1), one cylinder was fit to the diaphysis of the metatarsal, and a second was fit to the distal metatarsal articular surfaces. The axis of the second cylinder defined the anatomical medial-lateral axis in CYL, and a series of vector cross products between the cylinder axes produced the orthogonal coordinate system CYL, which is essentially aligned with the metatarsal head and shaft. CYL fitting was performed by a single operator in geomagic software (Geomagic Studio, 3D Systems, Rock Hill, SC). All proximal phalange coordinate systems were defined using principal components of the phalanges. Given the varying, oblate shape of the two sesamoids, both of their coordinate system orientation definitions were simply copied from the proximal phalange and moved to their respective volumetric centroids.

Six foot and ankle tendons responsible for great toe and ankle articulation and stabilization (Achilles, extensor digitorum longus, extensor hallucis longus, peroneus longus, peroneus brevis, and tibialis posterior) were dissected by an orthopedic surgeon, and sutured with surgical monofilament using Krackow stitches to interface via loops with braided aluminum tendon pull cables. Tibia and fibula shafts were exposed 60 mm proximal to the ankle mortise, denuded of tissues, and cast in resin cylinders. A custom loading rig allowing for four tibial segments kinematic degrees-of-freedom (ab/adduction and medial-lateral translation were fixed) was assembled around the WBCT scanner in a manner that did not interfere with scan quality. Given the lack of in vivo muscle force data for these activities, tendon excursions and loadings were extrapolated from previous cadaveric experiments mimicking foot function [39–41]. Static axial half-bodyweight donor-specific loads were applied through the tibial shafts with free weights. A computer-controlled actuator with a load cell mounted in series with the Achilles tendon induced the target ankle plantarflexion and consequent physiologic distal foot motion. All other tendons were nominally and statically loaded with 15 N.

Two motions were simulated: (1) an isolated toe lift (extension of the great toe) and (2) heel rise of the entire foot via Achilles tendon actuated ankle plantarflexion while the forefoot remains on the ground. For the toe lift, a custom foot plate made of acrylic and adjustable to specimen size was machined and instrumented with a potentiometer to control and record the precise amount of toe extension (Fig. 2). Toe extension targets ranged from 0 deg (neutral foot) to the maximum of the specimen (approximately 50 deg or 60 deg) in 10 deg increments. For each toe lift pose target, the extensor hallucis longus tendon was pulled to induce the motion, and the instrumented toe plate was used to set the exact rotation amount and lock the specimen into a supported position at the desired extension angle with the tendon clamped in tension. For the heel rise testing, specimens were positioned directly on the floor of the WBCT as living subjects would stand. By varying the Achilles load and allowing natural articulation of the joints and translation of the tibial shaft in the testing apparatus, seven scan positions equally spaced between a neutral foot flat pose and the maximum stable heel rise pose were achieved. Specimens were held for 30 s at each pose before scanning to allow tissues to reach a viscoelastic steady-state. If tissue creep or clamp slippages were detected during the postscan review, scans were recollected.

Volumetric rigid registration was performed between the specimen-specific template volumes and the experimental scans for each of the four bones to determine their positions and orientations across all poses. Registrations were performed in custom matlab (Mathworks, Natick, MA) software that optimally aligns the specimen- and bone-specific template volumes to the raw WBCT volumes using a cost function based on the mutual information between volumes [42,43]. This optimization metric proved to be robust to the noise present in our low-dose WBCT data. A validation of this registration process is presented in the Appendix.

Joint kinematics were derived using the transformations describing each template volume to the global CT frame and the local bone coordinate systems. Motions were expressed relative to the metatarsal CYL and PCA coordinate system, separately. The resulting matrix was decomposed into Euler angles in the following order: sagittal plane flexion–extension rotation, transverse plane ab/adduction rotation, and frontal plane inversion–eversion rotation.

Surface models of each bone were generated in Mimics and filtered with Laplacian smoothing to reduce segmentation noise. For the MTPJ1, medial MTSJ, and lateral MTSJ the weighted centroids of articular joint contact were estimated using the opposing subchondral bone surface distances at each reconstructed pose [44]. Centroid locations were combined over sequential poses to define contact paths of the proximal phalange and sesamoids on the metatarsal articular surface.

The finite helical axes (FHA) of motion were determined for the MTPJ1, medial MTSJ, and lateral MTSJ using the method described by Spoor and Veldpaus [45] and Woltring [46]. Using all permutations of pose pairings within each simulated activity, the mean FHA (MFHA) was determined using the symmetrical axis of rotation analysis approach described by Ehrig with a $w= sin(θ2)$ weighting term to minimize errors from small rotations [47,48]. This produced robust estimations of the mean functional axes of each joint for both simulated activities. The median locations of the instantaneous centers of rotation (ICR) in the sagittal plane were also calculated using the R package gmedian (v 1.2.5) [49] for MTPJ1 and each MTSJ [50,51].

Left foot specimen data were mirrored to appear as right, and all metatarsals were coregistered together iteratively in matlab. From these data, a mean metatarsal model was generated in the University of Utah's shapeworksstudio software [52,53]. In matlab, the affine transformations of each specimen's metatarsal to the volume representing the mean bone shape were determined and used to scale, mirror, and map the MFHA, and arthrokinematic data of each specimen onto the common metatarsal model for illustration and statistical comparison. The mean morphological flexion–extension axis was calculated across all PCA and CYL definitions. Angles of inclination of this flexion–extension axis relative to the MFHA were calculated in both the transverse plane (α) and frontal plane (β).

With all data mapped into a common reference frame centered in the mean metatarsal head, directional statistics [54] were utilized to quantify the spatial orientations of the morphological (PCA and CYL) and functional axes (MFHA) from each joint/activity (Fig. 3). Correlations were calculated between the MFHA vector spatial orientations and the first metatarsal medial-lateral axes for both PCA and CYL definitions (H1). To assess differences in the resultant kinematics between the PCA and CYL methods

(H2), a repeated measures ANOVA was performed across all poses with the axis definition type (PCA or CYL) as the within-subjects factor. Kinematic crosstalk was quantified by correlating the metatarsophalangeal flexion–extension data with the two remaining orthogonal rotational components. Analyses were performed in R statistical software using the base and directional (v4.8) packages [55,56]. This investigation tested the following hypotheses: (H1) MFHA of MTPJ1 and MTSJ will exhibit stronger spatial correlations to medial-lateral metatarsal axis definitions based on CYL compared to PCA, and (H2) that CYL-based coordinate systems will produce MTPJ1 kinematic data with less crosstalk and less intersubject variability compared PCA.

The biomechanical testing rig produced physiological foot poses for each specimen and activity. The Achilles tendon failed in the clamp for one specimen during the heel rise simulation, but enough tissue remained for reclamping and continued testing without noticeable effect. Repeatability of the CYL axis definitions was assessed via single-rater intraclass correlation for consistency of measurement (ICC = 0.985). Details of the bone tracking validation, and bias and precision data for each bone are presented in the Appendix.

Qualitatively, contact paths for each bone were consistent within specimens (low variability, data not shown) and grossly similar across activities (Fig. 4). MTPJ1 contact paths deviated slightly laterally in the heel rise activity, compared to the nearly vertical pattern observed in toe lift (Fig. 4). The lateral sesamoid path moved more medially in heel rise compared to toe lift. Qualitatively, the longest paths of contact were observed in the MTPJ1 compared to the sesamoids (Fig. 5). Heel rise activities induced longer contact paths than simple toe extension. In all conditions, the shortest contact paths were observed in the lateral sesamoid compartment.

Median ICRs for all bones were in the metatarsal head (Fig. 6). Sesamoid ICRs were located proximal and dorsal to the phalange ICRs, which were essentially centered in the distal metatarsal head. Qualitatively, proximal phalange ICRs were clustered in tighter distributions compared to both sesamoids for both activities. Compared to toe lift, the heel rise-activity sesamoid ICRs were shifted distally and dorsally (Fig. 6).

The orientations of the morphological axes and mean finite helical axes for each joint (metatarsophalangeal, lateral metatarsal-sesamoid, and medial metatarsal-sesamoid) were approximately normal to the sagittal plane (Fig. 7 and Table 1). For the PCA and CYL axes, there was no significant difference in the transverse plane orientation (mean difference 1.7 deg; p = 0.173), while the frontal plane orientation difference approached significance (mean difference 7.1 deg; p = 0.051). The medial sesamoid MFHA orientation was significantly different from the PCA and CYL in seven of eight comparisons across both planes and activities (Table 1). Qualitatively, in heel rise activities, the MFHA orientations were less variable in the anterior-posterior direction and spread more along the vertical (head-ground) direction (Fig. 7). The orientations of the MFHA were strongly correlated with both the PCA and CYL flexion–extension axes for both toe off and heel rise (Table 2).

Joint kinematics resulting from each of the two metatarsal coordinate system definitions (PCA and CYL) were similar in appearance for the primary flexion–extension motion across both activities, with a significant near-constant offset between the curves (p = 0.006, partial eta squared = 0.201) (Figs. 8(a) and 8(b)). Larger variability was observed between PCA and CYL kinematics in the motions about the two remaining orthogonal axes (Figs. 8(c), 8(d), 8(e), and 8(f)). For the ab/adduction axis, the two definitions were not different (p = 0.155, partial eta squared = 0.1), likely because the PCA method exhibited a large variance at each pose (Figs. 8(c) and 8(d)). For inversion–eversion (p = 0.041, partial eta squared = 0.117), the PCA group exhibited a larger variance and a bias of about 5 deg of inversion.

Correlation between flexion–extension and inversion–eversion rotations was moderate for CYL (R² = 0.30) and strong for PCA (R² = 0.51). Correlation between flexion–extension and ab/adduction was weak for both CYL (R² = –0.18) and PCA (R² = –0.10). Crosstalk between inversion–eversion and ab/adduction was moderate for CYL (R² = –0.37) and strong for PCA (R² = –0.54).

This study quantified the kinematics and arthrokinematics of the MTPJ1 and MTSJs during two relevant simulated foot motions and compared the orientations of the mean functional helical axes of motion (MFHA) to the primary flexion–extension axes in two methods (PCA and CYL) for embedding anatomical metatarsal coordinate systems. While these data are based on controlled cadaveric experiments, the resultant sesamoidal and phalangeal arthrokinematics for both simulated activities are qualitatively physiological (Figs. 4 and 5), and agree with the limited published data, as discussed below.

Our ex vivo kinematic data (Fig. 4) compare well with previous reports of proximal phalange motion and sesamoid motion in cadavers [8–10]. As also reported by Ahn et al., we observed the dorsal migration of the MTPJ1 contact region [10]. Shereff et al. reported that both sesamoids displaced 10–12 mm over an arc of motion of 111 deg, whereas we found 4–6 mm of sesamoid motion over 50–60 deg of arc (Fig. 5) [8]. Jamal et al. quantified ex vivo sesamoid motion through 70 deg of rotation; they reported comparable lateral sesamoid sagittal excursions (7.5 mm), but found more medial sesamoid excursion (14.7 mm) [9]. These discrepancies may be due to limitations in biofidelity of the applied loads to cadavers or in the methodology for quantifying and reporting the motion in the previous studies.

Analysis of various joints' ICRs during functional tasks has proven useful in quantifying unstable articulations and for designing optimal joint replacement hardware and installation procedures [31,57]. We found the median ICR locations of all bones agreed well with Shereff's report that they reside well within the confines of the metatarsal head in normal feet (Figs. 6 and 7) [8]. In hallux valgus and rigidus patients, Shereff reported that ICRs of the proximal phalange and sesamoids are located in more variable positions than those typically found in the “normal” cases, indicating how MTPJ1 pathology can disrupt joint kinematics. In a theoretical model, Durrant argued the ICRs of both the MTPJ1 and MTSJs should reside within the metatarsal head, and their distributions would be influenced by the cam shape of the distal head; specifically the flatter distal plantar surface [58]. The distributions (Fig. 7) of our proximal phalange ICR data fell roughly along the metatarsal long axes, while the sesamoid ICRs were distributed toward the dorsum. This may be indicative of the sesamoids having greater translation and sliding motion rather than pure rotation about the cam-shaped distal plantar surface compared to the proximal phalange's motion relative to the metatarsal. Further quantification of ICR locations and MFHA orientations may be clinically useful in diagnosing forefoot issues and in designing joint replacement strategies that respect the native relationships between these functional axes.

This study also compared the downstream effects of embedding two types of coordinate systems in the metatarsal. Our first hypothesis (H1), that spatial correlation to MFHA would be stronger for CYL than PCA was not supported (Table 2). Our second hypothesis (H2), that resultant kinematics from PCA would be more variable (Fig. 8) and have greater crosstalk compared to CYL, as indicated by stronger correlations between primary and secondary kinematic rotations with PCA, is currently supported by our data. While the PCA method uses the entire bone surface geometry and may account for gross morphological differences [59], natural variation in twisting of the metatarsal head relative to the diaphysis exists that may bias definitions of the medial-lateral axis at the MTPJ1. Kitashiro described a natural 11–13 deg mean bias toward eversion in first metatarsal torsion (range: –2 to 25 deg) [24]. Drapeau compared this torsion across primates and found a mean bias of 7.6±5.5 deg toward eversion for homo sapiens [37]. Both reports agree well with our observed PCA axis biases in the frontal plane (6.3±9.4 deg of eversion). Additionally, Mortier suggested that metatarsal torsional patterns and the orientation of the crista, which is responsible for guiding sesamoid motion, have direct connections to hallux valgus progression [38]. Furthermore, Brenner reported 7.99±3.71 deg of crista tilt relative to the metatarsal long axis, which may account for the larger transverse plane variances in sesamoid MFHA orientations in our study [36]. These studies suggest that more data characterizing the relationships between form and function in the MTPJ1 and MTSJs are required.

While great care was taken to simulate physiological motions in minimally dissected specimens, there are many limitations to our methodologies. Repeated trials were not collected for each specimen, the same tendon load scaling was applied to all specimens regardless of the sex or size of the donor, and the intrinsic muscles of the foot were not loaded in order to preserve specimen integrity against excessive dissection. Differences in kinematics observed between coordinate system types may be exacerbated or nullified in different tasks or with a larger pool of subjects, or by exploring alternative methods for defining the coordinate systems. We utilized an affine registration procedure to align all specimen data to a common mean metatarsal model, and while the vector orientation differences are not appreciably different when recalculated using rigid alignments (Appendix Table 3), the mean model to which each specimen was registered is a function of the study sample size and is assumed to represent a population mean. Given the presence of rotations in the other two planes, it might be an oversimplification to model the MTPJ1 joint as a hinge normal to the sagittal plane. While the robust symmetrical axis of rotation analysis approach was utilized to characterize MFHA, measurement noise and small rotations contribute to errors and the resultant axes are specific to the prescribed tasks and ranges of motion simulated [47,48]. Therefore, the efficacy of these methods must be tested during in vivo functional tasks like gait and stair navigation and in cohorts with atypical and highly varying bony morphologies.

In conclusion, our study outlines a protocol for quantitatively assessing the efficacy and utility of anatomically derived joint coordinate systems by comparing the geometric parameters of the articulations (axes of best-fit) to dynamic kinematic data (helical axes of motion). Our findings demonstrated the benefit of utilizing joint coordinate systems based on fitting geometric primitives to the underlying anatomic metatarsal bone geometry and articulations over the simpler method using the principal components of the entire bone. While embedding coordinate systems based on PCA is computationally straightforward to implement, the downstream effects of this decision on joint kinematic crosstalk errors must be considered. The relationships between the geometric and kinematic axes need to be further investigated in living-subjects performing multiple activities, and with various clinical diagnoses of MTPJ1 and sesamoid-related issues. Despite the challenges associated with imaging these small bones, their function and motion are essential to normal toe biomechanics, and more quantitative data on their motion during relevant tasks may inform and advance foot treatment and understanding.

• U.S. Department of Veterans Affairs (Grant Nos. RX002357 and RX002970; Funder ID: 10.13039/100000738).

• University of Washington Medical School Scholarship of Discovery Program.

CYL =

cylinder fitting-based coordinate system

ICR =

instantaneous center of rotation

MFHA =

mean finite helical axis

MTPJ1 =

first metatarsal phalangeal joint

MTSJ =

metatarsosesamoidal joint

PCA =

principal component analysis-based coordinate system

WBCT =

weightbearing computed tomography

α =

transverse plane inclination angle

β =

frontal plane inclination angle

We utilized an affine registration of each specimen's metatarsal model to the mean bone shape of all specimens. Affine registration allows for nonuniform scaling and shear in addition to rotations and translations. This morphs the specimen model, particularly the geometric relationships between the distal head of the metatarsal and its body toward the mean shape. It may also have the potential to artificially distort the differences between the various vector orientations examined in this study. We decomposed the affine transformation matrices used for each specimen into the rigid (translations and rotations) and nonrigid (shear and scale) matrices in matlab (qr.m function) and recomposed the transformations with only the rigid components. The vector orientations for the MFHA, and the CYL and PCA axes were remapped to the mean metatarsal model using these rigid transformations, and the angular deviations between each are reported in Table 3.

The accuracy and precision of reconstructing the first metatarsal, proximal phalanx, and medial and lateral sesamoid bones with our volumetric registration algorithm were assessed via a cadaveric validation study. Upon conclusion of the study, a specimen that best represented the median size and motion across all tested specimens was selected for further imaging in the WBCT scanner. A custom set of fiducial registration clusters fabricated consisting of carbon fiber shafts and Delrin spheres. Carbon fiber material was chosen for its lightweight and minimally attenuating properties; essential to avoiding artifacts during scanning or bending moments on the small sesamoid bones. The Delrin spheres (5 and 8 mm diameter) were mounted in triad configurations, and the base of each carbon shaft was rigidly embedded into each predrilled bone (Fig. 9). Plantar fat and skin were dissected to allow uninhibited movement of the carbon shafts when the bones were actuated via tendon pulling. WBCT scanning was performed with the great toe flexed and extended through its range of motion in five equally spaced positions. For each scan, the poses of the four bones of interest were reconstructed three separate times using our custom software utilizing volumetric tracking with a mutual information cost function [43]. The resultant poses were compared to the “gold standard” positions determined by segmenting the fiducial sphere clusters in mimics software (Materialize, Belgium) and determining the appropriate transformations. The translational and rotational bias and precision errors for each joint are listed in Table 4. Our bone registration algorithm could accurately track the translations of all the bones, and accurately determined the rotations of the metatarsal and phalanx. The rotational orientations of the sesamoids bones had larger errors, which was expected given their smaller bone volumes and oblate shapes. The limitations for this validation include the sample size (N = 1) and the inability to directly validate the heel rise activity, as our design for the sesamoid marker clusters could be constrained by the ground.