Multi-level multi-domain statistical shape model of the subtalar, talonavicular, and calcaneocuboid joints

Traditionally, two-dimensional conventional radiographs have been the primary tool to measure the complex morphology of the foot and ankle. However, the subtalar, talonavicular, and calcaneocuboid joints are challenging to assess due to their bone morphology and locations within the ankle. Weightbearing computed tomography is a novel high-resolution volumetric imaging mechanism that allows detailed generation of 3D bone reconstructions. This study aimed to develop a multi-domain statistical shape model to assess morphologic and alignment variation of the subtalar, talonavicular, and calcaneocuboid joints across an asymptomatic population and calculate 3D joint measurements in a consistent weightbearing position. Specific joint measurements included joint space distance, congruence, and coverage. Noteworthy anatomical variation predominantly included the talus and calcaneus, specifically an inverse relationship regarding talar dome heightening and calcaneal shortening. While there was minimal navicular and cuboid shape variation, there were alignment variations within these joints; the most notable is the rotational aspect about the anterior-posterior axis. This study also found that multi-domain modeling may be able to predict joint space distance measurements within a population. Additionally, variation across a population of these four bones may be driven far more by morphology than by alignment variation based on all three joint measurements. These data are beneficial in furthering our understanding of joint-level morphology and alignment variants to guide advancements in ankle joint pathological care and operative treatments.

1 Introduction

Historically, the complex nature of foot and ankle joint morphology has primarily been analyzed individually from two-dimensional (2D) measurements on conventional radiographs (Lopez-Ben, 2015; Krähenbühl et al., 2017). These methods fail to illustrate the complexity of the foot and ankle joints and their spatial relationships. The subtalar, talonavicular, and calcaneocuboid joints are specifically challenging to assess radiographically due to their intricate morphologies and locations within the ankle (Ebraheim et al., 1999; Willauer et al., 2014; Krähenbühl et al., 2017; Bernasconi et al., 2021). Yet, accurately visualizing these three joints is crucial to comprehending the compensatory joint mechanics as well as treatment for multiple pathologies, such as osteoarthritis and progressive collapsing flatfoot deformity (PCFD) (Sammarco, 2004; Bruening et al., 2018; Welte et al., 2021). However, advancements in volumetric imaging, including computed tomography (CT), have made it possible and practical to generate high-resolution three-dimensional (3D) reconstructions of bones throughout the foot and ankle (Hayes et al., 2006; Barg et al., 2018). While these imaging modalities are typically performed in a non-weightbearing position, weightbearing cone-beam CT (WBCT) technology allows for the analysis of joint relationships in a natural and consistent position with the presence of load (Colin et al., 2014; Burssens et al., 2016; Krähenbühl et al., 2016). A better understanding of hind- and midfoot morphometrics could assist in diagnosing and treating multiple joint diseases. However, it can be challenging to quantitatively extract morphological metrics from 3D surface reconstructions while maintaining anatomical relevance.

Statistical shape modeling (SSM) is a population-based mathematical approach to objectively quantify these morphological metrics (Davies et al., 2002; Styner et al., 2003; Cates et al., 2007; Datar et al., 2009; Datar et al., 2011; Goparaju et al., 2022). Using SSM, a statistical model can be created to compare mean bone shape morphology and identify anatomical modes of variation. Previous ankle SSM studies have been limited to a single bone and could not, therefore, evaluate multi-domain joint relationships (Melinska et al., 2015; Melinska et al., 2017; Tümer et al., 2019a; Tümer et al., 2019b; Gabrielli et al., 2020; Krähenbühl et al., 2020; Liu et al., 2020; Schmutz et al., 2021; Arbabi et al., 2022; Peiffer et al., 2022; Vafaeian et al., 2022). A multi-domain technique can be implemented to capture morphological and alignment changes for multiple bones throughout a population. Additionally, multi-level analyses allow for separating morphology and alignment to identify their individual contributions to joint relationships.

Using a multi-domain SSM approach, individual joint-level 3D morphometrics can be calculated to predict variations within the joints based on the original alignment from WBCT scans. These joint metrics include coverage area, joint space distance, and joint congruence. The coverage area is a mathematically derived calculation of the joint’s articulating region. Coverage area is also used to quantify morphology variation, such as osteoarthritis development (Schaefer et al., 2012), and alignment variation, such as joint subluxation (Louie et al., 2014). Calculated from within the joint’s coverage area, the joint space distance is calculated from the Euclidean distance across the joint, and the congruence index rates how well the two articular surfaces match one another (Ateshian et al., 1992; Lenz et al., 2021). Similar to coverage area variations, joint space distance can help clinically indicate degenerative diseases like osteoarthritis (Day et al., 2020) and pathologies with varying alignments, such as PCFD (Bernasconi et al., 2021). These analyses allow for a holistic understanding of bone relationships in the foot and ankle that can guide implant design and development while providing insight into the involvement of multiple joints within various pathological diseases.

This study aims to characterize asymptomatic joint-level morphology and alignment differences throughout statistically significant modes of variation and the joint-level measurements within the subtalar, talonavicular, and calcaneocuboid joints. Additionally, this study aims to determine which, if any, multi-level multi-domain SSM modes of variation that seek to separate morphology from joint alignment can be used to predict joint-level measurements across a population. The validity of the multi-domain SSM approach was assessed by comparing the multi-domain statistical models to models run via a single-domain SSM approach. We hypothesized that the first mode of variation could predict joint coverage, distance, and congruence for the population. We also hypothesized that multi-domain SSM does not affect the mean shape compared to single-domain SSM. To accomplish these aims and test these hypotheses, we have presented a new computation approach to evaluate multi-domain influences of morphology and alignment in the case of four bones within the ankle and statistically evaluated the downstream joint measurements.

2 Materials and methods

Twenty-seven asymptomatic participants (age: 50.0 ± 7.3 years; height: 169.4 ± 6.4 cm; BMI: 25.3 ± 3.8 kg/m²; seven males) were enrolled with IRB approval from an ethics committee (Kantonsspital Baselland, Switzerland; University of Utah). Individuals between 40 and 70 years of age without a history of ankle injury or surgery were considered. Before the study, a clinical and radiographic assessment was performed to exclude participants with a planovalgus or cavovarus deformity.

Each participant underwent a unilateral weightbearing CT (WBCT) scan (Planmed Verity, Planmed Oy, Helsinki, Finland; 0.4 mm isotropic pixel resolution). The WBCT scans were segmented, decimated, and smoothed to generate 3D surface models of the talus, calcaneus, navicular, and cuboid (Amira, v6.0.1, Visage Imaging, San Diego, CA, United States).

2.1 Statistical shape modeling

For each of the twenty-seven individuals, the talus, calcaneus, navicular, and cuboid were used to create two types of statistical models utilizing an opensource SSM software (ShapeWorks v6.2.1, University of Utah; www.shapeworks.sci.utah.edu). This constructed Particle-based Shape Models (PSM) that automatically placed a dense set of corresponding landmarks on the given set of shapes using an entropy-based optimization scheme (Cates et al., 2017). The SSM approaches included a single-domain SSM (SD SSM) of only the individual bones and a multi-domain SSM (MD SSM) comprised of all four bones in their anatomical alignment (Figure 1). Before optimizing the shape model, any bone models from a left limb were reflected to the right. Using a global iterative closest point algorithm, the input bones for each of the statistical models were aligned. The individual bones were aligned separately for the single-domain model, and each of the four bones were aligned together for the multi-domain model. Aligning the bone models together maintained the individual anatomical alignments from the WBCTs (Wilm, 2022).

FIGURE 1

FIGURE 1. Workflow terminology for the model analysis performed with the simplified terminology bolded. The multi-domain statistical shape model principal component analysis (MD SSM) is a four bone model maintaining anatomical joint relationships analyzing morphology and alignment. The multi-level principal component analysis in alignment shape space (alignment variant) is a four bone model maintaining anatomical joint relationships solely analyzing alignment. The multi-level principal component analysis in morphology shape space (morphology variant) is a four bone model maintaining anatomical joint relationships solely analyzing morphology. The single-domain statistical shape model mean validation (SD SSM) is four single bone models to validate morphology variation to the MD SSM.

The shape model for the four bones is built together by a multi-domain shape modeling approach (Cates et al., 2008), where point correspondences for all the surfaces are optimized in the full joint space of the object complex. A total particle count of 1,024 for the talus, 2,048 for the calcaneus, 512 for the navicular, and 512 for the cuboid was used for both the single and multi-domain statistical models. These correspondence particles were used to define mean shapes and statistically quantify shape differences between bones within the different statistical models. A Procrustes analysis was not applied for these statistical models to remove scale as a factor. Patient specific joint measurements include size variations across the population. In order to compare our population joint level measurements with SSM models, we chose not to use Procrustes for equivalent statistical comparisons. A principal component analysis (PCA) was used to simplify the data to a smaller set of linearly uncorrelated components, or modes of variation. Using a parallel analysis algorithm, statistically significant modes of variation were calculated (p < 0.05) (Horn, 1965; Ledesma and Valero-Mora, 2007). For each significant mode of variation, the surface distances between the mean and first and second standard deviation (SD) shapes were calculated and visualized via CloudCompare (v2.11. alpha, www.cloudcompare.org). The shape model was analyzed using two different approaches. The first is the multi-domain approach (MD SSM), where PCA analysis is done on the entire multi-bone joint structure treated as a single bone complex reflecting the modes of variation of the entire multi-object complex. The other analysis approach is the multi-level technique, where PCA analysis is done separately on the shape of each individual bone and on the alignment of the entire multi-bone joint structure. This multi-level approach allowed for separately analyzing the alignment shape space model (alignment variant) and the morphology shape space model (morphology variant) (Figure 1). The same significant modes of variation were used for subsequent calculations. For the morphology variant, surface distances between the mean shape and ±1 SD shapes were calculated and visualized via CloudCompare, similar to the MD SSM approach. And for the alignment variant, the mean shape and ±1 SD shapes for each mode were overlaid on each other to visualize those variations.

2.2 Joint coverage, distance, and congruence index calculations

The correspondence particles from each statistical model were used to automatically calculate and compare the joint-level measurements with an available toolbox (Peterson, 2022). Joint coverage was calculated automatically using normal vectors from the faces of the individual bone models. For a given bone of a joint, if the normal vectors from that bone intersected with a face of the opposing bone, they were considered within the coverage region for that side of the joint. The surface area of the coverage region was then calculated by summing the surface area of the faces that comprised the coverage region. Subject-specific correspondence particles within the coverage region would then be used for that bone’s side of the joint analysis. These identified correspondence particles were paired with their nearest neighboring surface mesh node. It was from this paired node that the joint space distance and congruence index calculations were made. This was necessary as calculating the congruence index requires the mean and Gaussian curvatures of the two opposing surfaces. The joint space distance at a correspondence particle was calculated as the Euclidean distance between its paired node and the nearest opposing surface node. The congruence index for that correspondence particle was then calculated between these two nodes. Congruence index, first described by Ateshian et al., is a rating of how well two surfaces match one another, with a congruence index of 0 mm^-1 being rated as perfectly congruent (Ateshian et al., 1992; Lenz et al., 2021). Both results were mapped to that specific correspondence particle. For visualization, correspondence particles with a distance value greater than 6 mm were colored white for both joint space distance and congruence index. Subsequent statistics were performed on all correspondence particles within the coverage regions

2.3 Population vs. PCA modal comparison

The mean and ± 1 SD for coverage across the twenty-seven participants for all four joints were automatically calculated. These average values will be referred to as population calculations throughout this study. Due to the nature of correspondence particles, each distance value and congruence index were paired with identical correspondence particles across the entire population. Thus, the mean and ± 1 SD for distance and congruence can be calculated at each correspondence particle

These calculations for each measurement across the population were statistically compared to the mean and ±1 SD for each significant mode of variation using an unpaired t-test based on the means, SDs, and population size. This statistical comparison was also used to compare the alignment and morphology of multi-level multi-domains to the population. A p-value less than 0.05 indicates that the two groups are significantly different from one another.

3 Results

3.1 Statistical shape model

Statistical analysis is performed using Principal Component Analysis (PCA) for the MD SSM, morphology variant and alignment variant, where the mean shape and modes of variation are computed based on the optimized MD SSM (Figure 1). This study observed three statistically significant PCA modes of variation, accounting for 74.8% of the overall shape variation. The three modes represented 63.8%, 6.3%, and 4.7% of the variation. Additionally, when comparing the SD SSM of each bone to the MD SSM, there were negligible differences (<0.1 mm) in mean shape using the previously mentioned surface distance calculations in CloudCompare. Computationally, this statistical shape model took about 40 min to run optimization and the joint measurement analysis took an average of about 3 min to run per patient for the largest joint.

3.1.1 Multi-domain statistical shape model analysis (MD SSM)

Substantial morphological and configurational variations were observed while maintaining joint articular relationships within the MD SSM approach. The first PCA mode of variation highlighted size variation primarily, with overall growth and shrinkage of all four bones simultaneously (Figure 2). The second PCA mode is multifaceted in the presented variation. The primary regions of change occur in the inverse relationship between the talus and calcaneus. Concerning the talus, as the talar dome heightens, the posterior process diminishes. Concerning the calcaneus, as the calcaneus lengthens, the posterior facet’s slope decreases. Analyzing the two bones simultaneously, as the talar dome heightens and the posterior process diminishes, the calcaneus shortens, and the posterior facet’s slope increases (Figure 2). Additionally, as the talar dome heightens, the navicular and cuboid move inferiorly with minimal rotation about the anterior-posterior axis. The third PCA mode primarily varies with the anteromedial facet moving anterior to posterior throughout the SDs. There is minor talar dome variation but little else of note within this general multi-domain approach (Figure 2).

FIGURE 2

FIGURE 2. MD SSM modes of variation one to three showing both morphological and alignment significant variations at ± 1 and 2 standard deviations from the mean shape. Black arrows highlight key anatomical feature variations. Red regions are larger than the mean shape and blue regions are smaller than the mean shape with a scale in millimeters (mm).

3.1.2 Multi-level analysis in alignment shape space (alignment variant)

When focusing on the alignment variant, notable anatomical alignment variations were observed. The first PCA mode of variation shows an overall outward and inward movement between the bones, effectively increasing and decreasing joint space distance (Figure 3). The second PCA mode primarily highlights the superior and inferior motion of the four bones. Specifically, as the talus moves superiorly, the calcaneus, cuboid, and navicular move inferiorly (Figure 3). And the third PCA mode identifies the medial and lateral movement of the talus and calcaneus; as the talus moves medially, the calcaneus moves laterally. Additionally, as the talus moves medially, the navicular rotates superior and lateral, and the cuboid rotates inferior and medial (Figure 3).

FIGURE 3

FIGURE 3. Alignment variant modes of variation one to three showing significant alignment variations at ± 2 standard deviations from the mean shape. Black arrows highlight key anatomical feature variations. Dark purple regions are +2 standard deviations and pink regions are -2 standard deviations from the mean shape in gray.

3.1.3 Multi-level analysis in morphology shape space (morphology variant)

With just observing the morphological variant, the first PCA mode of variation shows each bone growing and shrinking individually (Figure 4). The second PCA mode of variation primarily has the lengthening of the calcaneus with a decreasing posterior facet slope. There are still slight talar dome changes and slight navicular and cuboid changes, but they are not as prominent as the general multi-domain approach (Figure 4). The third PCA mode shows a similar anterior/posterior anteromedial facet variation; however, it illustrates more of a rotational component. The anteromedial facet’s slope changes throughout the SDs from a steep slope to a more flattened slope (Figure 4).

FIGURE 4

FIGURE 4. Morphology variant modes of variation one to three showing significant morphology variations at ± 1 and 2 standard deviations from the mean shape. Black arrows highlight key anatomical feature variations. Red regions are larger than the mean shape and blue regions are smaller than the mean shape with a scale in millimeters (mm).

3.2 Joint measurements

3.2.1 Coverage

3.2.1.1 Subtalar joint

The coverage area for the subtalar joint was calculated for the entire subtalar joint and the posterior and anteromedial facets. The coverage area across the population for the entire joint averaged 1,519.34 ± 193.22 mm² for the talus and 1,435.99 ± 171.36 mm² for the calcaneus. None of the modes of variation within the MD SSM, morphology variant, or alignment variant were significantly different from the population. Moreover, the talus had a consistently larger coverage area than the calcaneus (Table 1).

TABLE 1

TABLE 1. Average coverage surface area ± standard deviation (SD) for the subtalar joint (split between the posterior and anteromedial facet), talonavicular joint, and calcaneocuboid joint on both involved surfaces. p-values are reported from each mode and variant to the population, and p-values reported between each bone within each joint. Bolded p-values are statistically significant to the population and between each bone, respectively.

The coverage area across the population for the posterior facet averaged 805.39 ± 112.42 mm² for the talus and 862.60 ± 127.85 mm² for the calcaneus. Similarly to the entire subtalar joint, none of the modes of variation within the MD SSM, morphology variant, or alignment variant differed significantly from the population. However, unlike the entire subtalar joint, the calcaneus had a consistently larger coverage area than the talus (Table 1).

The coverage area across the population for the anteromedial facet averaged 365.28 ± 74.76 mm² for the talus and 384.33 ± 77.94 mm² for the calcaneus. While none of the modes of variation within the MD SSM, morphology variant, or alignment variant for the talus were significantly different from the population, most of them were different when comparing the calcaneus to the population. Additionally, the talus coverage area is generally larger compared to the calcaneus coverage area, even though that is not true for the population (Table 1).

3.2.1.2 Talonavicular joint

The coverage area across the population for the talonavicular joint averaged 471.22 ± 67.08 mm² for the talus and 541.36 ± 55.15 mm² for the navicular. Mode one for the MD SSM, morphology variant and alignment variant was the only mode that was significantly similar to the population for the talar coverage area, each with p-values greater than 0.05. None of the p-values showed significant differences for the MD SSM, morphology variant or alignment variant compared to the population (Table 1). Additionally, the navicular had a significantly larger coverage area when compared to the talus in every analysis, including the population.

3.2.1.3 Calcaneocuboid joint

The coverage area across the population for the calcaneocuboid joint average 393.92 ± 55.07 mm² for the calcaneus and 406.46 ± 53.86 mm² for the cuboid. For the calcaneus, only the first mode of the morphology variant was significantly similar compared to the population. For the cuboid, only the first modes of the MD SSM and morphology variant were significantly similar compared to the population. When comparing the coverage area of the calcaneus to the cuboid, most were significantly different (Table 1). Still, the population, MD SSM mode one, and all modes in the morphology variant were not significantly different between the bones. However, the cuboid had a larger coverage area when compared to the calcaneus for every analysis, including the population.

3.2.2 Joint space distance

3.2.2.1 Subtalar joint

The population across the entire subtalar joint averaged a joint space distance of 3.33 ± 2.06 mm, with a narrower joint space in the posterior and lateral regions of the posterior facet and the medial region of the anteromedial facet (Figure 5). The average joint space distance throughout the first three modes of the MD SSM, alignment variant and morphology variant ranged from 3.40 ± 2.10 mm to 3.54 ± 2.12 mm (Table 2). There were minor variations across the joints throughout the second and third modes of variation. And while the first mode of variation had large changes across the first SDs, none of the subtalar distance values had a p-value lower than 0.7109 compared to the population (Table 2).

FIGURE 5

FIGURE 5. Average joint space distance with ±1 standard deviation for the subtalar articular region. Results are visualized for the population, MD SSM modes one to three, alignment variant modes one to three and morphology variant modes one to three via correspondence particles. Results are reported in millimeters (mm). Values larger than 6 mm in joint space distance are colored white.

TABLE 2

TABLE 2. Average joint space distance ±standard deviation (SD) for the subtalar joint, talonavicular joint, and calcaneocuboid join. p-values are reported from each mode and variant to the population. Bolded p-values are statistically significant to the population.

3.2.2.2 Talonavicular joint

The average talonavicular joint space distance for the population was 1.32 ± 0.43 mm, with a slight widening of the medial and central aspects of the joint (Figure 6). The average joint space distance throughout the first three modes of the MD SSM, alignment variant and morphology variant ranged from 1.27 ± 0.83 mm to 1.35 ± 0.88 mm (Table 2). While the first mode expectedly had larger changes across the first SD, the second mode also had substantial variation in both morphology and alignment variants. However, they all had a similar pattern of widening distance towards the central aspect of the joint, and none had a p-value lower than 0.7822 compared to the population (Table 2).

FIGURE 6

FIGURE 6. Average joint space distance with ±1 standard deviation for the talonavicular articular region. Results are visualized for the population, MD SSM modes one to three, alignment variant modes one to three and morphology variant modes one to three via correspondence particles. Results are reported in millimeters (mm). Values larger than 6 mm in joint space distance are colored white.

3.2.2.3 Calcaneocuboid joint

The population across the calcaneocuboid joint averaged a joint space distance of 1.67 ± 0.62 mm with a narrowing in the inferolateral region of the joint (Figure 7). The average joint space distance throughout the first three modes of the MD SSM, alignment variant and morphology variant ranged from 1.64 ± 0.37 mm to 1.73 ± 1.13 mm (Table 2). Similar to the previous joints, the first mode has the widest variation across the SD, with modes two and three showing a similar pattern as the mean. Additionally, none of the calcaneocuboid distance values had a p-value lower than 0.8098 compared to the population (Table 2).

FIGURE 7

FIGURE 7. Average joint space distance with ±1 standard deviation for the calcaneocuboid articular region. Results are visualized for the population, MD SSM modes one to three, alignment variant modes one to three and morphology variant modes one to three via correspondence particles. Results are reported in millimeters (mm). Values larger than 6 mm in joint space distance are colored white.

3.2.3 Joint congruence index

3.2.3.1 Subtalar joint

The population across the entire subtalar joint had an average congruence index of 0.24 ± 0.16 mm⁻¹. With a worsening congruence index on the outer edge of the joint (Figure 8). The average joint congruence index throughout the first three modes of the MD SSM, alignment variant and morphology variant ranged from 0.16 ± 0.14 mm⁻¹ to 0.20 ± 0.18 mm⁻¹ (Table 2). All modes and variants had similar congruence patterns and none had a p-value lower than 0.0559 compared to the population (Table 3). Overall, the population congruence was consistently higher than the MD SSM, morphology variant and alignment variant.

FIGURE 8

FIGURE 8. Average joint congruence indices with ±1 standard deviation for the subtalar articular region. Results are visualized for the population, MD SSM modes one to three, alignment variant modes one to three and morphology variant modes one to three via correspondence particles. Results are reported in inverse millimeters (mm⁻¹). Values larger than 6 mm in joint space distance are colored white.

TABLE 3

TABLE 3. Average congruence index ±standard deviation (SD) for the subtalar joint, talonavicular joint, and calcaneocuboid join. p-values are reported from each mode and variant to the population. Bolded p-values are statistically significant to the population.

3.2.3.2 Talonavicular joint

The average talonavicular joint congruence index across the population was 0.23 ± 0.12 mm^−1, with subjectively consistent congruence across the joint (Figure 9). The average joint congruence index throughout the first three modes of the MD SSM, alignment variant and morphology variant ranged from 0.16 ± 0.06 mm⁻¹ to 0.22 ± 0.07 mm⁻¹ (Table 3). While all the modes and variants had similar congruence patterns, the only variant with a p-value larger than 0.05 was the first mode in the morphology variant (0.7099). All other congruence values were significantly different in a more congruent fashion compared to the population. Overall, the population congruence was consistently higher than the MD SSM, morphology variant and alignment variant.

FIGURE 9

FIGURE 9. Average joint congruence indices with ±1 standard deviation for the talonavicular articular region. Results are visualized for the population, MD SSM modes one to three, alignment variant modes one to three and morphology variant modes one to three via correspondence particles. Results are reported in inverse millimeters (mm⁻¹). Values larger than 6 mm in joint space distance are colored white.

3.2.3.3 Calcaneocuboid joint

The average calcaneocuboid joint congruence index across the population was 0.30 ± 0.18 mm⁻¹ with worsening congruence towards the edge of the joint (Figure 10). The average joint congruence index throughout the first three modes of the MD SSM, alignment variant and morphology variant ranged from 0.21 ± 0.12 mm⁻¹ to 0.22 ± 0.15 mm⁻¹ (Table 3). Only the first and second MD SSM modes and the first alignment variant mode had a p-value larger than 0.05. All other congruence values were significantly different in a more congruent fashion compared to the population. Overall, the population congruence was consistently higher than the MD SSM, morphology variant and alignment variant.

FIGURE 10

FIGURE 10. Average joint congruence indices with ±1 standard deviation for the calcaneocuboid articular region. Results are visualized for the population, MD SSM modes one to three, alignment variant modes one to three and morphology variant modes one to three via correspondence particles. Results are reported in inverse millimeters (mm⁻¹). Values larger than 6 mm in joint space distance are colored white.

4 Discussion

The study’s primary aim was to use multi-domain SSM to characterize asymptomatic joint-level morphology and alignment variations throughout a population and determine joint-level measurements within the subtalar, talonavicular, and calcaneocuboid joints. Using those population-based joint-level measurements, we also aimed to determine if SSM can predict joint-level measurements. Moreover, we aimed to determine the validity of the multi-domain mean shapes compared to the single-domain mean shapes. The most relevant findings include: I) joint space distance is the only joint-level measurement statistically similar across all modes of variation for the MD SSM, morphology variant and alignment variant; II) the first mode in the morphology variant is statistically similar to the population across all joint-level measurements except the coverage area on the calcaneal anteromedial facet; III) MD SSM does preserve the mean shape when compared to the SD SSM mean shape.

Our MD SSM shape variations finding agree with the limited studies performed on the talus, calcaneus, navicular and cuboid. We found that the primary anatomical morphological variations included the talar trochlea height, calcaneal lengthening, and minimal variation throughout the navicular and cuboid (Melinska et al., 2015; Melinska et al., 2017; Tümer et al., 2019a; Krähenbühl et al., 2020; Lenz et al., 2021). However, we also found the relationship between the talar trochlea and calcaneal length to be interrelated. Additionally, the novel multi-level approach allowed for the separation of alignment from morphology. This approach provides insight of the talonavicular and calcaneocuboid alignment variations. Primary alignment variation occurred between the talus, navicular and cuboid. Specifically, when the talus moves inferiorly, the navicular and cuboid move superiorly and vice versa. But when the talus moves medially or laterally, the navicular and cuboid rotate about the anterior-posterior axis. Multiple studies have attempted to analyze and understand the role these midfoot joints have on human locomotion but have reached a wide range of conclusions. Some studies have found the talonavicular and calcaneocuboid joints to rotate throughout gait creating what is known as a locking and unlocking mechanism based on the axes of the transverse tarsal joint in one plane (Elftman, 1960; Suckel et al., 2008; Sarrafian, 2011); yet others have determined the motion within these joints did not significantly increase during gait (Blackwood et al., 2005). Moreover, other studies have determined that the plantar fascia has more of a role via the windlass mechanism in producing movement during gait (Hicks, 1954; Welte et al., 2021) and that there was not any observed midtarsal movement during gait (Bruening et al., 2018). There is still much to be understood about these different mechanisms and the role that the hind- and mid-foot joints have during gait, but multi-domain models may be a step towards understanding how morphology and alignment during static positioning influence these well-debated topics clinically.

Few studies have performed weightbearing joint space distance calculations of these three joints. One study that has calculated the asymptomatic subtalar joint distance for the joint as a whole reported 3.29 ± 0.87 mm compared to our measurement of 3.33 ± 2.06 mm, and when compared determined to not be statistically different from one another (p = 0.95) (Dibbern et al., 2021). While a small number of studies have calculated talonavicular and calcaneocuboid joint space distance, the most comparable to this study also used WBCT scans (Bernasconi et al., 2021). Bernasconi et al. reported average talonavicular and calcaneocuboid distance measurements of 1.14 ± 0.49 mm and 1.67 ± 1.26 mm, respectively. Compared to our measurements of 1.32 ± 0.43 mm for the talonavicular and 1.67 ± 0.62 mm for the calcaneocuboid, these measurements were statistically similar to one another (talonavicular: p = 0.28 and calcaneocuboid: p = 1.0). As previously mentioned, the first relevant finding of this study was that joint space distance was the only joint-level measurement statistically similar between the population and all modes of variation for the MD SSM, morphology variant and alignment variant. This finding may indicate that MD SSM, morphology variant and alignment variant models have the potential to predict joint distance measurements within a population.

When separating alignment (alignment variant) from morphology (morphology variant) in the first mode of variation, the relationship between size and joint space can be determined. As seen in the SSM figures, as the models get larger towards -2 SDs (Figure 4), they also become further apart (Figure 3). This relationship can be further seen in the joint space distance figures (Figures 5–7) with the morphology indicating that joint space is narrowing as the bone models get larger, but the alignment indicating that they are simultaneously moving further apart. This pattern indicates that individuals with larger bones subsequently have greater joint space. While joint space has been previously correlated to sex, height and weight measurements (Goker et al., 2009), further correlation to bone size and shape in the foot and ankle may be concluded from this approach. Additionally, the second relevant finding referring to the first mode of variation in the morphology variant being generally statistically similar to the population gives potential relevance to another predictive conclusion. With the first mode of the morphology variant producing the most statistically similar joint measurement results to the population, it may indicate that this specific mode can be used to predict all three measurements within a population. It also may indicate that variation across a population of these four bones is driven far more by morphology than by alignment variation. To further support this postulate, identifying patient populations with morphological variations, such as those with asymmetric ankle osteoarthritis, those suffering from PCFD, or after fracture reduction, could help better understand the role morphology plays in disease and deformity progression.

This study is not without its limitations. Soft tissues, specifically ligaments, tendons, and articular cartilage, were not included in this study. All joint-level measurements were calculated on subchondral bone surfaces, which do not consider the effects of these soft tissues when concerning alignment and joint distance. Second, WBCT scans are dependent on the native CT resolution and therefore maximum precision is 0.4 mm when considering CT segmentations for joint measurements. However, with the error spectrum being ±1 pixel, it is still less than 1 mm, therefore we feel confident in our ability to report clinically meaningful joint measurements. Third, static imaging does not capture the articulation behavior of these joints that imaging during human locomotion could capture. Further analyses on joint-level measurements during dynamic activities should be evaluated to further understand the true relationship these joints have with one another. Fourth, this study only includes bone models from healthy asymptomatic individuals. While this was intentionally done to determine healthy variations, incorporating pathological bone models in future studies could give better context of the morphological variations found in this healthy model. Further studies relating these findings to individuals with hindfoot ankle diseases will provide more relevance to clinical actions that should be taken.

In conclusion, joint-level morphology and alignment variations can be further understood using a multi-level multi-domain SSM. This study found joint space distance measurements across all modes of variation and the first mode of the morphology variation across all primary measurements to be statistically similar to the population. This may indicate using SSM to predict joint-level measurements in specific variants is a valid approach. Additionally, multi-domain SSM does preserve mean shape compared to a single-domain SSM mean shape and can confidently be used to further multi-domain SSM studies. Further studies include expanding multi-domain modeling to subsequent joints in the foot and ankle for asymptomatic and symptomatic populations. Additionally, performing multi-domain joint-level measurements during dynamic activities is necessary to fully understand the complex relationship between these joints in the foot and ankle. Future optimization schemes could be developed where the morphology and alignment variants are coherently reflected in the particle position updates of the optimized correspondence model.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by University of Utah Institutional Review Board. The patients/participants provided their written informed consent to participate in this study.

Author contributions

AP, NiK, CS, AB, and AL conceived and designed the experiment, NiK provided image data, AP, RL and AL generated surface models for all bones, NaK and SE developed statistical shape model software and multi-level analysis and provided technical oversight and troubleshoot on the use of the different modelling and analysis methods for the proposed study, AP conducted the statistical shape model development and analysis, AP, RL and AL conducted the joint analysis for coverage, distance, and congruence, and AP analyzed the results. All authors reviewed the manuscript.

Funding

The National Institutes of Health supported this work under grant numbers NIAMS-K01AR080221, NIBIB-U24EB029011, NIAMS-R01AR076120, and NIBIB-R01EB016701.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Author disclaimer

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

References

Arbabi, S., Seevinck, P., Weinans, H., Jong, P. A., Sturkenboom, J., Hamersvelt, R. W., et al. (2022). Statistical shape model of the talus bone morphology: A comparison between impinged and nonimpinged ankles. J. Orthop. Res. doi:10.1002/jor.25328

CrossRef Full Text | Google Scholar

Ateshian, G. A., Rosenwasser, M. P., and Mow, V. C. (1992). Curvature characteristics and congruence of the thumb carpometacarpal joint: Differences between female and male joints. J. Biomech. 25, 591–607. doi:10.1016/0021-9290(92)90102-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Barg, A., Bailey, T., Richter, M., De Cesar Netto, C., Lintz, F., Burssens, A., et al. (2018). Weightbearing computed tomography of the foot and ankle: Emerging technology topical Review. Foot Ankle Int. 39, 376–386. doi:10.1177/1071100717740330

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernasconi, A., De Cesar Netto, C., Siegler, S., Jepsen, M., and Lintz, F. (2021). Weightbearing CT assessment of foot and ankle joints in Pes Planovalgus using distance mapping. Foot Ankle Surg. 28, 775–784. doi:10.1016/j.fas.2021.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Blackwood, C. B., Yuen, T. J., Sangeorzan, B. J., and Ledoux, W. R. (2005). The midtarsal joint locking mechanism. Foot Ankle Int. 26, 1074–1080. doi:10.1177/107110070502601213

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruening, D. A., Pohl, M. B., Takahashi, K. Z., and Barrios, J. A. (2018). Midtarsal locking, the windlass mechanism, and running strike pattern: A kinematic and kinetic assessment. J. Biomech. 73, 185–191. doi:10.1016/j.jbiomech.2018.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Burssens, A., Peeters, J., Buedts, K., Victor, J., and Vandeputte, G. (2016). Measuring hindfoot alignment in weight bearing CT: A novel clinical relevant measurement method. Foot Ankle Surg. (N. Y). 22, 233–238. doi:10.1016/j.fas.2015.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Cates, J. E., Elhabian, S. Y., and Whitaker, R. T. (2017). ShapeWorks: Particle-Based shape correspondence and visualization software. Stat. Shape Deformation Analysis 2017, 257–298. doi:10.1016/B978-0-12-810493-4.00012-2

CrossRef Full Text | Google Scholar

Cates, J., Fletcher, P. T., Styner, M., Hazlett, H. C., and Whitaker, R. (2008). Particle-based shape analysis of multi-object complexes. Med. Image Comput. Comput. Assist. Interv.International Conf. Med. Image Comput. Computer-Assisted Intervention 11, 477–485. doi:10.1007/978-3-540-85988-8_57

PubMed Abstract | CrossRef Full Text | Google Scholar

Cates, J., Fletcher, P. T., Styner, M., Shenton, M., and Whitaker, R. (2007). Shape modeling and analysis with entropy-based particle systems. Inf. Process. Med. Imaging 20, 333–345. doi:10.1007/978-3-540-73273-0_28

PubMed Abstract | CrossRef Full Text | Google Scholar

Colin, F., Horn Lang, T., Zwicky, L., Hintermann, B., and Knupp, M. (2014). Subtalar joint configuration on weightbearing CT scan. Foot Ankle Int. 35, 1057–1062. doi:10.1177/1071100714540890

PubMed Abstract | CrossRef Full Text | Google Scholar

Datar, M., Cates, J., Fletcher, P. T., Gouttard, S., Gerig, G., and Whitaker, R. (2009). Particle based shape regression of open surfaces with applications to developmental neuroimaging. Med. Image Comput. Comput. Assist. Interv. 12, 167–174. doi:10.1007/978-3-642-04271-3_21

PubMed Abstract | CrossRef Full Text | Google Scholar

Datar, M., Gur, Y., Paniagua, B., Styner, M., and Whitaker, R. (2011). Geometric correspondence for ensembles of nonregular shapes. Med. Image Comput. Comput. Assist. Interv. 14, 368–375. doi:10.1007/978-3-642-23629-7_45

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, R. H., Twining, C. J., Cootes, T. F., Waterton, J. C., and Taylor, C. J. (2002). A minimum description length approach to statistical shape modeling. IEEE Trans. Med. Imaging 21, 525–537. doi:10.1109/tmi.2002.1009388

PubMed Abstract | CrossRef Full Text | Google Scholar

Day, M. A., Ho, M., Dibbern, K., Rao, K., An, Q., Anderson, D. D., et al. (2020). Correlation of 3D joint space width from weightbearing CT with outcomes after intra-articular calcaneal fracture. Foot Ankle Int. 41, 1106–1116. doi:10.1177/1071100720933891

PubMed Abstract | CrossRef Full Text | Google Scholar

Dibbern, K. N., Li, S., Vivtcharenko, V., Auch, E., Lintz, F., Ellis, S. J., et al. (2021). Three-dimensional distance and coverage maps in the assessment of peritalar subluxation in progressive collapsing foot deformity. Foot Ankle Int. 42, 757–767. doi:10.1177/1071100720983227

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebraheim, N. A., Haman, S. P., Lu, J., and Padanilam, T. G. (1999). Radiographic evaluation of the calcaneocuboid joint: A cadaver study. Foot Ankle Int. 20, 178–181. doi:10.1177/107110079902000307

PubMed Abstract | CrossRef Full Text | Google Scholar

Elftman, H. (1960). The transverse tarsal joint and its control. Clin. Orthop. 16, 41–46.

PubMed Abstract | Google Scholar

Gabrielli, A. S., Gale, T., Hogan, M., and Anderst, W. (2020). Bilateral symmetry, sex differences, and primary shape factors in ankle and hindfoot bone morphology. Foot Ankle Orthop. 5, 247301142090879. doi:10.1177/2473011420908796

PubMed Abstract | CrossRef Full Text | Google Scholar

Goker, B., Gonen, E., Demirag, M. D., and Block, J. A. (2009). Quantification of the radiographic joint space width of the ankle. Clin. Orthop. Relat. Res. 467, 2083–2089. doi:10.1007/s11999-009-0832-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Goparaju, A., Iyer, K., Bône, A., Hu, N., Henninger, H. B., Anderson, A. E., et al. (2022). Benchmarking off-the-shelf statistical shape modeling tools in clinical applications. Med. Image Anal. 76, 102271. doi:10.1016/j.media.2021.102271

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayes, A., Tochigi, Y., and Saltzman, C. L. (2006). Ankle morphometry on 3D-CT images. Iowa Orthop. J. 26, 1–4.

PubMed Abstract | Google Scholar

Hicks, J. H. (1954). The mechanics of the foot. II. The plantar aponeurosis and the arch. J. Anat. 88, 25–30.

PubMed Abstract | Google Scholar

Horn, J. L. (1965). A rationale and test for the number of factors in factor Analysis. Psychometrika 30, 179–185. doi:10.1007/bf02289447

PubMed Abstract | CrossRef Full Text | Google Scholar

Krähenbühl, N., Horn-Lang, T., Hintermann, B., and Knupp, M. (2017). The subtalar joint: A complex mechanism. EFORT Open Rev. 2, 309–316. doi:10.1302/2058-5241.2.160050

PubMed Abstract | CrossRef Full Text | Google Scholar

Krähenbühl, N., Lenz, A. L., Lisonbee, R. J., Peterson, A. C., Atkins, P. R., Hintermann, B., et al. (2020). Morphologic analysis of the subtalar joint using statistical shape modeling. J. Orthop. Res. 38, 2625–2633. doi:10.1002/jor.24831

PubMed Abstract | CrossRef Full Text | Google Scholar

Krähenbühl, N., Tschuck, M., Bolliger, L., Hintermann, B., and Knupp, M. (2016). Orientation of the subtalar joint: Measurement and reliability using weightbearing CT scans. Foot Ankle Int. 37, 109–114. doi:10.1177/1071100715600823

PubMed Abstract | CrossRef Full Text | Google Scholar

Ledesma, R. D., and Valero-Mora, P. (2007). Determining the number of factors to retain in EFA: An easy-to-use computer program for carrying out parallel analysis. Pract. Assess. Res. Eval. 12, 234. doi:10.7275/wjnc-nm63

CrossRef Full Text | Google Scholar

Lenz, A. L., Krähenbühl, N., Peterson, A. C., Lisonbee, R. J., Hintermann, B., Saltzman, C. L., et al. (2021). Statistical shape modeling of the talocrural joint using a hybrid multi-articulation joint approach. Sci. Rep. 11, 7314. doi:10.1038/s41598-021-86567-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Jomha, N. M., Adeeb, S., El-Rich, M., and Westover, L. (2020). Investigation of the average shape and principal variations of the human talus bone using statistic shape model. Front. Bioeng. Biotechnol. 8, 656. doi:10.3389/fbioe.2020.00656

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopez-Ben, R. (2015). Imaging of the subtalar joint. Foot Ankle Clin. 20, 223–241. doi:10.1016/j.fcl.2015.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Louie, P. K., Sangeorzan, B. J., Fassbind, M. J., and Ledoux, W. R. (2014). Talonavicular joint coverage and bone morphology between different foot types. J. Orthop. Res. 32, 958–966. doi:10.1002/jor.22612

PubMed Abstract | CrossRef Full Text | Google Scholar

Melinska, A. U., Romaszkiewicz, P., Wagel, J., Antosik, B., Sasiadek, M., and Iskander, D. R. (2017). Statistical shape models of cuboid, navicular and talus bones. J. Foot Ankle Res. 10, 6. doi:10.1186/s13047-016-0178-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Melinska, A. U., Romaszkiewicz, P., Wagel, J., Sasiadek, M., and Iskander, D. R. (2015). Statistical, morphometric, anatomical shape model (atlas) of calcaneus. PLoS One 10, e0134603. doi:10.1371/journal.pone.0134603

PubMed Abstract | CrossRef Full Text | Google Scholar

Peiffer, M., Burssens, A., De Mits, S., Heintz, T., Van Waeyenberge, M., Buedts, K., et al. (2022). Statistical shape model-based tibiofibular assessment of syndesmotic ankle lesions using weight-bearing CT. J. Orthop. Res. doi:10.1002/jor.25318

CrossRef Full Text | Google Scholar

Peterson, A. C. (2022). Static joint measurement analysis. [Online] Available: https://github.com/Lenz-Lab/StaticJointMeasureAnalysis.

Google Scholar

Sammarco, V. J. (2004). The talonavicular and calcaneocuboid joints: Anatomy, biomechanics, and clinical management of the transverse tarsal joint. Foot Ankle Clin. 9, 127–145. doi:10.1016/s1083-7515(03)00152-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarrafian, S. (2011). Sarrafian's anatomy of the foot and ankle. Philadelphia, Pennsylvania, United States: Lippincott Williams & Wilkins.

Google Scholar

Schaefer, K. L., Sangeorzan, B. J., Fassbind, M. J., and Ledoux, W. R. (2012). The comparative morphology of idiopathic ankle osteoarthritis. J. Bone Jt. Surgery-American Volume 94, e181. doi:10.2106/jbjs.l.00063

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmutz, B., Lüthi, M., Schmutz-Leong, Y. K., Shulman, R., and Platt, S. (2021). Morphological analysis of Gissane’s angle utilising a statistical shape model of the calcaneus. Arch. Orthop. Trauma Surg. 141, 937–945. doi:10.1007/s00402-020-03566-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Styner, M., Gerig, G., Lieberman, J., Jones, D., and Weinberger, D. (2003). Statistical shape analysis of neuroanatomical structures based on medial models. Med. Image Anal. 7, 207–220. doi:10.1016/s1361-8415(02)00110-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Suckel, A., Muller, O., Langenstein, P., Herberts, T., Reize, P., and Wulker, N. (2008). Chopart's joint load during gait. Gait Posture 27, 216–222. doi:10.1016/j.gaitpost.2007.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Tümer, N., Arbabi, V., Gielis, W. P., De Jong, P. A., Weinans, H., Tuijthof, G. J. M., et al. (2019a). Three-dimensional analysis of shape variations and symmetry of the fibula, tibia, calcaneus and talus. J. Anat. 234, 132–144. doi:10.1111/joa.12900

PubMed Abstract | CrossRef Full Text | Google Scholar

Tümer, N., Vuurberg, G., Blankevoort, L., Kerkhoffs, G. M. M. J., Tuijthof, G. J. M., and Zadpoor, A. A. (2019b). Typical shape differences in the subtalar joint bones between subjects with chronic ankle instability and controls. J. Orthop. Res. 37, 1892–1902. doi:10.1002/jor.24336

PubMed Abstract | CrossRef Full Text | Google Scholar

Vafaeian, B., Riahi, H. T., Amoushahi, H., Jomha, N. M., and Adeeb, S. (2022). A feature-based statistical shape model for geometric analysis of the human talus and development of universal talar prostheses. J. Anat. 240, 305–322. doi:10.1111/joa.13552

PubMed Abstract | CrossRef Full Text | Google Scholar

Welte, L., Kelly, L. A., Kessler, S. E., Lieberman, D. E., D'andrea, S. E., Lichtwark, G. A., et al. (2021). The extensibility of the plantar fascia influences the windlass mechanism during human running. Proc. R. Soc. B 288, 20202095. doi:10.1098/rspb.2020.2095

PubMed Abstract | CrossRef Full Text | Google Scholar

Willauer, P., Sangeorzan, B. J., Whittaker, E. C., Shofer, J. B., and Ledoux, W. R. (2014). The sensitivity of standard radiographic foot measures to misalignment. Foot Ankle Int. 35, 1334–1340. doi:10.1177/1071100714549188

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilm, J. (2022). Iterative closest point. New York: MATLAB Central File Exchange.

Google Scholar