Elastic Modulus and Its Relation to Apparent Mineral Density in Juvenile Equine Bones of the Lower Limb

Limb fractures in horses often result in euthanasia due to limitations of internal fixation related to body mass and anatomy. Up to 80% of racehorse fatalities are caused by a fracture [1], and this number has not improved since the mid-1970s [2]. The majority of fatal musculoskeletal injuries in the lower limb of racing horses occur in the third metacarpus (MC3) and proximal phalanx (P1) [3,4] and are the result of chronic fatigue [2]. Epidemiological studies have linked several factors to increased fracture risk including racetrack surface, injury history, and sex [1], and significant effort has been placed toward addressing environmental risk factors. While some progress has been made, the goal of preventing essentially all fractures has yet to be realized.

Bone is a functionally adaptive material that responds to its local mechanical environment [5]. Exercise in young horses, while the skeleton is primed for adaptation, has been shown to increase P1 diaphyseal bone mineral content and bone area [6], suggesting an opportunity to direct bone modeling in such a way to reduce fracture risk later in life. Computational models can be used to noninvasively predict the mechanical loading environment of bone in vivo [7,8] and therefore provide a means for evaluating the effect of different exercise regimens preclinically rather than adopting a trial and error approach. Critical to these predictions are accurate material properties, such as the Young's modulus, which can be related empirically to computed tomography (CT) based apparent mineral density [9,10].

Several CT density–modulus relationships exist for horse bone [10–12]. However, age of the samples ranged from 3 to 14 years, representing adult bone since equine skeletal maturity occurs at approximately 2 years of age. The most widely used density–modulus relationship for equine bone was developed by Les et al. [10] for longitudinal MC3 samples, but the mean sample age was 6.7 years. The nonlinear nature of the density-modulus relationship makes it difficult to determine whether one may extrapolate existing functions to young bone. Moreover, structural and compositional differences between juvenile (immature) and adult (mature) bone also motivate the need to consider a different density–modulus relationship for foals than for adult horses.

Like other large mammals, young horses initially have plexiform cortical bone, which contains woven bone, that eventually converts to Haversian bone as they mature [13]. Chappard and colleagues also described juvenile trabecular bone as “plexiform” but this has not been reported elsewhere [14]. In ovine trabecular bone, mature bone has increased bone volume fraction and apparent ash density, and decreased collagen content when compared to immature bone [15,16]. Similarly, elastic modulus, ultimate stress, and ultimate strain are known to be different in juvenile bone compared to adult bone [15,17].

In human trabecular bone, it has been established that density–modulus relationships vary by anatomical location [18] and are anisotropic [9], but whether this is true in equine trabecular bone is not known. For example, the growth plates in the distal MC3 and proximal P1 close at the same time [19], however whether or not the bones mineralize at the same rate is not known. Therefore, we hypothesize that juvenile equine bone may require a different density–modulus relationship than those currently reported. Thus, the objective of this study was to develop a density–modulus relationship for juvenile equine bone and evaluate the sensitivity of the density–modulus relationship to anatomical location and loading direction.

Intact bones (n = 18) were collected from young horses euthanized for reasons unrelated to this study (Table 1). Distal limbs were collected within 4 h of euthanasia and frozen at −20 °C. Prior to subsequent steps, distal limbs were cleaned of soft tissue, disarticulated, and the MC3 and P1 were wrapped in PBS soaked gauze and stored in sealed plastic bags at −20 °C.

The bones and four mineral density phantoms (range: 25–750 mg HA/cm³, CIRS, Norfolk, VA) were scanned in a clinical CT scanner (LightSpeed16, GE HealthCare, Chicago, IL) with the same protocols used for live horses (nominal voxel resolution = 0.875 × 0.875 × 0.625 mm, 120 kVp, 200 mA). To avoid artifacts (overestimates of apparent CT density) associated with scanning excised cores [20], we developed a method to identify the bone cores virtually within the intact image volume. First, the intact bones were imaged with overlapping micro-CT (μCT) scans (nominal isotropic resolution = 144 μm, 90 kVp, 177 μA, Rigaku CTLab GX130, Tokyo, Japan) acquired along the bone length and merged. The intact clinical and μCT scans were aligned in three-dimensional space to share a coordinate system origin and slice plane (Fig. 1(a)). The bone sections that remained following core extraction were μCT scanned and aligned to the intact CT data to locate each core within the μCT dataset (Fig. 1(b)). Virtual location of the cores in the μCT dataset was then transformed to the location in the clinical CT dataset (Fig. 1(c), Amira 3D v2022.2, Thermo Fisher Scientific, Waltham, MA). Custom Matlab code (v2022b, MathWorks, Natick, MA) was then used to separate core CT stacks from the entire bone CT stack via a masking process.

In order to maximize the number of cores that could be extracted from each bone, the intact μCT was used to assess the trabecular structure and randomly assign sections of the bone (perpendicular to the long axis) to either longitudinal or transverse cores. Each bone section was cut using a water-irrigated diamond band saw (Model C-40 Tall, Gryphon, Sylmar, CA) while the bone was frozen. Trabecular cores were removed in the longitudinal (MC3 n = 85; P1 n = 49) and dorsal-palmar transverse (MC3 n = 54; P1 n = 36) directions using a water-irrigated diamond sintered coring bit (internal diameter = 5 mm, His Glassworks, Asheville, NC) mounted on a drill press. The ends of the cores were trimmed and ground perpendicular to the long axis using sandpaper wetted with PBS (grit: 220, 500, 800). Bone marrow was left intact [21], and the cores were fixed in custom Delrin endcaps using two-part epoxy (endcap diameter = 19 mm, endcap length $≈$ 10 mm, exposed length = 12.36 ± 1.33 mm, total embedded length = 4.83 ± 1.87 mm, Fig. 2(a)). Custom jigs were used during the potting process to ensure the bone sample remained perpendicular to the plane of the endcaps. Between core machining and embedding in endcaps the samples were wrapped in PBS soaked gauze and stored individually in 5 mL Eppendorf tubes at −20 °C. After embedding in endcaps, samples were wrapped in PBS-soaked gauze and refrigerated at approximately 2 °C for 18 h before testing.

Cores were returned to room temperature and thoroughly hydrated with PBS. Compression tests were performed on a tabletop test frame (Instron 5967, Norwood, MA) using a fixed lower platen and a self-aligning, spherical seated upper platen. Embedded cores were preloaded to 5N, preconditioned for five cycles by loading to 0.001 strain at 0.01 strain/s, and then loaded at 0.01 strain/s until failure [10,21]. Force and crosshead displacement were recorded at 100 Hz. Crosshead displacement measures both system compliance and sample displacement [21]. Therefore, we determined an adjusted crosshead displacement measure that accounted for machine compliance (using the direct technique [22]) and compliance of the endcap material.

The diameter and the exposed length of each core was measured using digital calipers. Stress was calculated by dividing force by each sample's cross-sectional area (diameter = 4.97 ± 0.04 mm). Strain was calculated by dividing adjusted crosshead displacement by effective gauge length (exposed length + 1/2 total length embedded in endcaps [23]) of each sample. Young's modulus was calculated as the slope of the linear regression of all data between two points on the elastic portion (approximately 0.003–0.018 ε) of the stress–strain curve.

Normality of Young's modulus within each anatomical location and orientation was evaluated using a Shapiro–Wilks test. Modulus in the MC3 longitudinal and transverse directions and P1 longitudinal direction were not normally distributed; therefore, distributions were compared between all groups using a Mann–Whitney–Wilcoxon test. To assess the effect of anatomical location and orientation on the density-modulus relationship, additional models were created with those variables as fixed effects, allowing for interaction with CT density, and compared to a model without the variable in question via likelihood ratio test to obtain a p-value. Linear mixed effects models do not have an R² in the traditional sense, therefore the method defined by Nakagawa and Schielzeth was used to calculate a marginal R² (⁠$Rm2$⁠) that represents the variance explained by the fixed factor (CT density) [25]. Conditional R² (⁠$Rc2$⁠), which represents the variance in modulus explained by CT density and random effect of subject combined, was also calculated for each relationship [25]. All analyses were performed in R (v4.2.1) and the lme4 package was used to perform the linear mixed effects analysis.

All data are presented as median ± median absolute deviation. Of the Young's modulus data shown in Fig. 3, only the P1 transverse data were normally distributed. The longitudinal modulus was 134% higher than the transverse modulus in the MC3 (Fig. 3(a)) and 90% higher in the P1 (Fig. 3(b)); these distributions differed significantly (p <0.001 for both). The transverse modulus of MC3 samples was lower (260 ± 230 MPa) than the P1 (351 ± 202 MPa, Fig. 3), although these distributions did not differ significantly (p = 0.12). There was no significant difference between longitudinal samples from the MC3 and P1 (609 ± 297 MPa and 667 ± 405 MPa, respectively).

Within the MC3 samples, orientation significantly affected model predictions of modulus (p <0.001), indicating that the longitudinal and transverse directions should have separate equations. When all P1 data were pooled, orientation significantly affected the model (p <0.001), again indicating models should be orientation specific. When data were pooled for each direction and anatomical location was included as a fixed effect, location significantly affected the model (p <0.001 in the longitudinal direction, p = 0.02 in the transverse direction). Together, these model results indicate that each anatomical location and orientation requires a different density–modulus relationship (Fig. 4).

Overall, the transverse modulus was better predicted than longitudinal modulus in both bones, with CT density predicting 86% of the variability in the P1 (⁠$EP1−T$⁠) and 77% of the variability in the MC3 (⁠$EMC3−T$⁠, Fig. 4). Variability in the longitudinal modulus was similarly predicted in both bones (⁠$Rm2$=0.62 for $EMC3−L$⁠, $Rm2$=0.66 for $EP1−L$⁠). Root mean squared percent error (RMSPE) was calculated for each of the density–modulus relationships to assess the magnitude of error in relation to actual values. In the MC3 (Fig. 4(a)), RMSPE was 35% in the longitudinal direction and 45% in the transverse direction. In the P1 (Fig. 4(b)), RMSPE was 32% in the longitudinal direction and 23% in the transverse direction.

Apparent CT density was converted to ash density using Eq. (2) in order to compare our juvenile MC3 longitudinal samples (E_young) to adult MC3 longitudinal samples (E_adult [12], Fig. 5). Root mean squared error (RMSE) was used to compare the fit of each model to the juvenile data points. The RMSE was 347 MPa using E_adult, compared to 191 MPa using E_young.

Using a robust sample size, we have developed the first density-modulus relationships for the juvenile equine MC3 and P1. To evaluate whether indeed this relationship is different from those derived from older bone, we compared our juvenile longitudinal MC3 data to that reported by Les et al. [10], and found that use of the adult relationship resulted in an almost 80% increase in error of modulus prediction. Although several equations relating hydroxyapatite derived CT density to ash density exist [24,26–28], we found they all led to an overprediction of modulus when using the adult density-modulus relationship [10] on juvenile samples. Modulus values for adult MC3 trabecular bone in compression range from 2.09 to 3.65 GPa [11,12,29], while median modulus for our juvenile MC3 longitudinal samples was 609 ± 297 MPa. It should be noted that over 80% of the samples tested by Les et al. were cortical bone (although they reported no difference in the density–modulus relationship between cortical and trabecular bone [10]), while our samples were trabecular. Bone volume fraction in the distal MC3 condyles reportedly increases from approximately 32% in 1–2 month old horses to approximately 60% in horses greater than 6 years old [30]. As apparent CT density is a combined measure of bone volume fraction and tissue density, it is expected that apparent CT density changes with maturation.

While differences in intrinsic properties between immature and mature bone are the most likely explanations for why density–modulus relationships differ with age, there are also methodological differences in the measurement of Young's modulus that are worth mentioning. We tested our samples using endcaps while Les et al. [10] tested samples in compression using platens directly in contact with the bone sample, which has since been shown to result in an underestimation of modulus by 20–40% [23] in trabecular bone samples. The underestimation of modulus would likely cause an even larger disparity between the density-modulus developed here for juvenile bone than that reported for older bone but this remains to be confirmed.

We found that density–modulus relationships in juvenile equine trabecular bone vary depending on anatomical location (MC3 versus P1) and loading direction (longitudinal versus transverse), which is consistent with human data [9,18]. The impact of not accounting for differences between bones when predicting modulus from CT data is notable. Applying the MC3 density–modulus relationships to the P1 data results in an RMSPE of 49% in the longitudinal direction and 31% in the transverse direction, which are both higher than the percent error obtained with using the P1 density–modulus relationships (longitudinal: 32%, transverse: 23%).

There are several reasons why the density–modulus relationships between bones may be different. Biomechanically, the two bones are likely under different types of loads with the MC3 in combined bending and compression due to its long slender structure and distribution of cortical material properties [31]. In contrast, the cuboidal P1 would be more resistant to bending. Variation in surface strain modes between the distal MC3 (compression) and the P1 (shear) have been reported [32,33].

We also found that the density–modulus relationships in the transverse direction had stronger predictions of modulus than the longitudinal direction (Fig. 4). The density–modulus relationship of the P1 in the transverse direction (⁠$EP1−T$⁠) had the highest $Rm2$ (0.86) and lowest percent error (RMSPE = 23%) of all relationships investigated. The density–modulus relationship of the MC3 in the transverse direction (⁠$EMC3−T$⁠) also had a high $Rm2$ (0.77) but had the highest percent error (RMSPE = 45%) of all relationships, which was driven by increased variability of the modulus data between a CT density of 0.2–0.5 g HA/cm³. For example, at a density of approximately 0.28 g HA/cm³ in the MC3 (Fig. 4(a)), the transverse modulus ranges from approximately 200–800 MPa (blue triangles).

The sensitivity of the strength of modulus predictions to orientation may be related to the nature of the microstructure along each direction. Longitudinal cores tend to have more varied microstructure along the length of the sample (Fig. 6(a)) when compared to the more compact microstructure evident in transverse cores (Fig. 6(b)). However, these microstructural differences may be unique to juvenile animals, as Keyak et al. reported similar R² values in density–modulus relationships for adult human proximal tibia bone in the longitudinal (0.84) and transverse directions (anterior–posterior: 0.72; medial-lateral: 0.84) [9]. Augat et al. found approximately equivalent coefficient of variation in modulus between the longitudinal and both transverse directions of adult human trabecular bone in the spine, calcaneus, proximal femur, and distal femur [34]. In addition to microstructure, mineral quantity, crystallinity, and carbonate substitutions change with subject age and are known to affect the elastic properties of bone [35,36]. Further work is needed to confirm whether mineralization rates are different between the MC3 and P1, as well as the influence of microstructure and tissue mineral density on the elastic modulus.

Although a driving motivator for this study was the ability to more accurately evaluate the mechanical environment during exercise in juvenile horses, there are other applications. Finite element models based on CT data have been used to assess fracture risk and represent an improved assessment of bone strength when compared to densitometric variables alone [37]. Finite element models can also be used as a presurgical evaluation tool to predict tissue response to certain fixation methods [38]. Use of the incorrect density–modulus relationship may result in incorrect assessments of bone strength or presurgical evaluation; therefore, our findings may have direct clinical implications.

There are several limitations of this study. Transverse bone samples were only machined in the dorsal-palmar (“anterior– posterior”) direction. The geometry of the MC3 and P1 led to challenges in excising cores in the medial-lateral direction, and limited availability of intact juvenile bones required we waste as little tissue as possible. The implications of this may be mitigated by the fact that trabecular bone has been described as a transversely isotropic structure [39], and modulus values in the anterior-posterior and medial-lateral directions often have a similar relationship with ash density [9] and bone volume fraction [40]. Despite the fact that we tested over 200 samples, we still encountered variability in the density–modulus data that leaves approximately 25% of the modulus variability unexplained in the case of the longitudinal MC3 (Fig. 4(a), $Rc2$ = 0.74 for $EMC3−L$⁠). Aside from the influence of microstructure on mechanical properties, some of this variability may be due to sample age, as we had donors ranging in age from approximately 0.5 week to 48 weeks. As well, 35% of samples in the MC3 and 60% of samples in the P1 (combining both orientations) were from subjects that were 4 weeks old (Table 1). Sample size in the current study does not allow us to develop statistically meaningful density-modulus relationships for each age of juvenile horses, and instead data were pooled to describe horses less than 1 year old. Nonetheless, all model predictions were significant and this work represents the first large scale mechanical testing study in juvenile equine bone.

Therefore, using rigorous imaging and experimental protocols, we have established orientation-specific density–modulus relationships for the juvenile MC3 and P1 bones. The incorporation of these data into computational models will allow for more accurate predictions of the mechanical response of young bone to loads and therefore the potential for bone adaptation.

Samples were collected by AMM during 2018–2022 from patients requiring euthanasia at the Large Animal Clinic of the Veterinary Teaching Hospital at UIUC. We are grateful for the technical support from Dr. Leilei Yin for his help with μCT scanning, and to the Visualization Laboratory at the Beckman Institute for Advanced Science and Technology for providing software used in this project. Funding was provided by USDA Hatch Funds (ILLU-888-964), Campus Research Board at UIUC, Grayson Jockey Club Research Foundation, and a graduate fellowship from the Beckman Institute for Advanced Science and Technology to SGM.

• Grayson-Jockey Club Research Foundation (Two year award granted in 2021; Funder ID: 10.13039/100001655).

• U.S. Department of Agriculture (Hatch Funds (ILLU-888-964); Funder ID: 10.13039/100000199).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.