Density-modulus relationships are necessary to develop ﬁnite element models of bones that may be used to evaluate local tissue response to di↵erent physical activities. It is unknown if juvenile equine trabecular bone may be described by the same density-modulus as adult equine bone, and how the density-modulus relationship varies with anatomical location and loading direction. To answer these questions, trabecular bone cores from the third metacarpal (MC3) and proximal phalanx (P1) bones of juvenile horses (age <1 yr) were machined in the longitudinal (n=134) and transverse (n=90) directions and mechanically tested in compression. Elastic modulus was related to apparent CT density of each sample using power law regressions. We found that density-modulus relationships for juvenile equine trabecular bone were signiﬁcantly di↵erent for each anatomical location (MC3 vs P1) and orientation (longitudinal vs transverse). Use of the incorrect density-modulus relationship resulted in increased root mean squared percent error of the modulus prediction by 8-17%. When our juvenile density-modulus relationship was compared to one of an equivalent location in adult horses, the adult relationship resulted in an 80% increase in error of the modulus prediction. Moving forward, more accurate models of young bone can be developed and used to evaluate potential exercise regimens designed to encourage bone adaptation. Keywords: juvenile, equine, density, modulus, third metacarpal (MC3), proximal phalanx (P1) 1. Introduction Limb fractures in horses often result in eu- thanasia due to limitations of internal ﬁxation related to body mass and anatomy. Up to 80% of racehorse fatalities are caused by a fracture 5 [1], and this number has not improved since the mid-1970s [2]. The majority of fatal musculoskele- tal injuries in the lower limb of racing horses oc- cur in the third metacarpus (MC3) and proximal phalanx (P1) [3, 4] and are the result of chronic 10 ⇤Corresponding author Email address: mkersh@illinois.edu (Mariana E. Kersh)fatigue [2]. Epidemiological studies have linked several factors to increased fracture risk includ- ing racetrack surface, injury history, and sex [1], and signiﬁcant e↵ort has been placed towards ad- dressing environmental risk factors. While some 15 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 20 primed for adaptation, has been shown to increase P1 diaphyseal bone mineral content and bone area [6], suggesting an opportunity to direct bone mod- eling in such a way to reduce fracture risk later in life. Computational models can be used to non- 25 invasively predict the mechanical loading environ- ment of bone in vivo [7, 8] and therefore provide a means for evaluating the e↵ect of di↵erent exer- cise regimens pre-clinically rather than adopting a trial and error approach. Critical to these pre- 30 dictions is accurate material properties, such as the Young’s modulus, which can be related em- pirically to computed tomography (CT) based ap- parent mineral density [9, 10]. Several CT density-modulus relationships exist 35 for horse bone [10, 11, 12]. However, the age of the samples ranged from 3-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 40 was developed by Les et al. [10] for longitudi- nal MC3 samples but the mean sample age was 6.7 years. The non-linear nature of the density- modulus relationship makes it dicult to deter- mine whether one may extrapolate existing func- 45 tions to young bone. Moreover, structural and compositional di↵erences between juvenile (im- mature) and adult (mature) bone also motivate the need to consider a di↵erent density-modulus relationship for foals than for adult horses. 50 Like other large mammals, young horses ini- tially have plexiform cortical bone, which contains woven bone, that eventually converts to Haver- sian bone as they mature [13]. Chappard and colleagues also described juvenile trabecular bone 55 as ‘plexiform’ but this has not been reported else- where [14]. In ovine trabecular bone, mature bone has increased bone volume fraction and appar- ent ash density, and decreased collagen content when compared to immature bone [15, 16]. Sim- 60 ilarly, elastic modulus, ultimate stress, and ulti- mate strain are known to be di↵erent in young bone compared to adult bone [15, 17]. In human trabecular bone it has been estab- lished that density-modulus relationships vary by 65 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 min- 70 eralize at the same rate is not known. There-fore, we hypothesize that juvenile equine bone may require a di↵erent density-modulus relation- ship than those currently reported. Thus, the objective of this study was to develop a density- 75 modulus relationship for juvenile equine bone and evaluate the sensitivity of the density-modulus re- lationship to anatomical location and loading di- rection. 2. Materials and Methods 80 2.1. Specimens Intact bones (n=18) were collected from young horses euthanized for reasons unrelated to this study (Table 1). Distal limbs were collected within 4 hours of euthanasia and frozen at 20C. 85 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 20C. Age (wk) 0.43 4 18 23 48 Bone MC3 P1 MC3 P1 MC3 P1 MC3 P1 MC3 P1 Intact Bones 11 37 11 11 11 Long. Cores 19 3 29 30 12 2 15 9 10 5 Trans. Cores 64 19 22 92 83 12 5 Table 1: Distribution of samples included in this study. Cores were removed in the longitudinal and transverse di- rections. 2.2. Sample Imaging 90 The bones and mineral density phantoms (range: 25-750 mg HA/cm3,C I R S )w e r es c a n n e d in a clinical CT scanner (LightSpeed16, GE Medi- cal Systems) with the same protocols used for live horses (nominal voxel resolution=0.875 x 0.875 95 x 0.625 mm, 120 kVp, 200 mA). To avoid arti- facts (overestimates of apparent CT density) as- sociated with scanning excised cores [20], we de- veloped a method to identify the bone cores vir- tually within the intact image volume. First, 100 the intact bones were imaged with overlapping microCT ( µCT) scans (nominal isotropic resolu- tion = 144 µm, 90 kVp, 177 µA, Rigaku CT- Lab GX130) acquired along the bone length and merged. The intact clinical and µCT scans were 105 aligned in 3D space to share a coordinate sys- tem origin and slice plane (Fig. 1A). The bone 2 Figure 1: Bone samples were collected from the MC3 and P1 bone of the equine forelimb (far left). (A) Clinical CT scan of a P1 (orange) aligned to whole bone µCT (dark blue) to share a coordinate system origin and slice plane. (B) Bone slabs scanned and aligned to whole bone µCT to locate cores virtually. (C) Bone core location identiﬁed within the clinical CT scan and density sampled within the core. sections that remained following core extraction were µCT scanned and aligned to the intact CT data to locate each core within the µCT data set 110 (Fig. 1B). Virtual location of the cores in the µCT data set were then transformed to the location in the clinical CT data set (Amira 2020.1, Fig. 1C). Custom Matlab code (v2021.b) was then used to separate core CT stacks from the entire bone CT 115 stack via a masking process. 2.3. Bone Core Preparation In order to maximize the number of cores that could be extracted from each bone, the intact µCT was used to assess the trabecular struc- 120 ture and randomly assign sections of the bone (perpendicular to the long axis) to either longi- tudinal or transverse cores. Each bone section was cut using a water-irrigated diamond band saw while the bone was frozen. Trabecular cores 125 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 di- amond sintered coring bit (internal diameter = 5mm) mounted on a drill press. The ends of the 130 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 ﬁxed in custom Delrin end- caps using 2-part epoxy (endcap diameter = 19 135 mm, endcap length ⇡10 mm, exposed length = 12.36 ±1.33 mm, total embedded length = 4.83 ± 1.87 mm). Custom jigs were used during the pot- ting process to ensure the bone sample remained perpendicular to the plane of the endcaps. Be- 140 tween core machining and embedding in endcaps Figure 2: (A) Schematic of sample embedded in Delrin endcaps (dimensions are not to scale). (B) Image of sample in the mechanical test setup. the samples were wrapped in PBS soaked gauze and stored individually in 5 mL Eppendorf tubes at -20C. After embedding in endcaps, samples were wrapped in PBS-soaked gauze and refrig- 145 erated at approximately 2Cf o r1 8h o u r sb e f o r e testing. 2.4. Compression Testing Cores were returned to room temperature and thoroughly hydrated with PBS. Compression 150 tests were performed on a tabletop test frame (In- stron 5967) using a ﬁxed lower platen and a self- aligning, spherical seated upper platen. Embed- ded cores were pre-loaded to 5N, pre-conditioned for 5 cycles by loading to 0.001 strain at 0.01 155 strain/s, and then loaded at 0.01 strain/s until failure. Force and crosshead displacement were 3 recorded after adjusting for machine compliance using the direct technique [22]. Measured dis- placement was also adjusted to account for end- 160 cap material compliance. Force and displacement data were sampled at 100 Hz. 2.5. Young’s Modulus Calculation The diameter and the exposed length of each core was measured using digital calipers. Stress 165 was calculated by dividing force by each sample’s cross-sectional area (diameter=4.97 ±0.04mm). Strain was calculated by dividing displacement by e↵ective gauge length (exposed length + 1/2 length embedded in endcaps [23]) of each sam- 170 ple. 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 ")o ft h es t r e s s - s t r a i nc u r v e . 2.6. Density measurement and density-modulus 175 relationship The clinical CT scan and phantoms were used to calculate average CT density ( ⇢CT, g HA/cm3) for each core. Modulus and CT density data were pooled by bone type (MC3, P1) and orientation 180 of core (longitudinal, transverse). In order to sat- isfy the assumptions of linearity, homoscedastic- ity, and normality of residuals, modulus and CT density were both log transformed. We used a linear mixed e↵ects regression between modulus 185 (dependent variable) and CT density (ﬁxed ef- fect), with subject included as a random e↵ect (random intercept). The slope and intercept were used to deﬁne the exponential and leading coe- cient terms, respectively, in the function relating 190 modulus to CT density. 2.7. Statistical Analysis Normality of Young’s modulus within each bone and anatomical location was evaluated using a Shapiro-Wilks test. Modulus in the MC3 lon- 195 gitudinal and transverse directions and P1 longi- tudinal direction were not normally distributed; therefore, distributions were compared between all groups using a Mann-Whitney-Wilcoxon test. To assess the e↵ect of anatomical location and 200loading direction on the density-modulus relation- ship, additional models were created with those variables as ﬁxed e↵ects, allowing for interaction with CT density, and compared to a model with- out the variable in question via likelihood ratio 205 test to obtain a p-value. Linear mixed e↵ects models do not have an R2in the traditional sense, therefore the method deﬁned by Nakagawa and Schielzeth was used to calculate a marginal R2 that represents the variance explained by the ﬁxed 210 factors (CT density) [24]. All analyses were per- formed in R (v4.2.1) and the lme4 package was used to perform the linear mixed e↵ects analysis. 3. Results Figure 3: Young’s modulus for the (A) MC3 and (B) P1 samples in the longitudinal (orange) and transverse (blue) directions. Distributions of mod- ulus were signiﬁcantly di↵erent between orientations within a bone, but not between bones. All data are presented as median ±median ab- 215 solute deviation. Of the Young’s modulus data 4 shown in Fig. 3, only the P1 transverse data was normally distributed. The longitudinal mod- ulus was 134% higher than the transverse mod- ulus in the MC3 (Fig. 3A) and 90% higher in 220 the P1 (Fig. 3B); these distributions di↵ered sig- niﬁcantly (p <0.001 for both). The transverse modulus of MC3 samples was lower (260 ±230 MPa) than the P1 (351 ±202 MPa, Fig. 3), al- though these distributions did not di↵er signiﬁ- 225 cantly (p=0.12). There was no signiﬁcant di↵er- ence between longitudinal samples from the MC3 and P1 (609 ±297 MPa and 667 ±405 MPa, respec- tively). Within the MC3 samples, orientation signif- 230 icantly a↵ected model predictions of modulus (p<0.001), indicating that the longitudinal and transverse directions should have separate equa- tions. When all P1 data were pooled, orientation signiﬁcantly a↵ected the model (p <0.001), again 235 indicating models should be orientation speciﬁc. When data were pooled for each direction and anatomical location was included as a ﬁxed e↵ect, location signiﬁcantly a↵ected the model (p <0.001 in the longitudinal direction, p=0.02 in the trans- 240 verse direction). Together, these model results indicate that each anatomical location and orien- tation requires a di↵erent density-modulus rela-tionship (Fig. 4). Overall, the transverse modulus was better pre- 245 dicted than longitudinal modulus in both bones, with CT density predicting 86% of the variability in the P1 and 77% of the variability in the MC3 (Fig. 4, blue lines). Variability in the longitudi- nal modulus was similarly predicted in both bones 250 (R2=0.62 in MC3, R2=0.66 in P1). 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. 4A), RMSPE was 35% 255 in the longitudinal direction and 45% in the trans- verse direction. In the P1 (Fig. 4B), RMSPE was 32% in the longitudinal direction and 23% in the transverse direction. 4. Discussion 260 Using a robust sample size, we have developed the ﬁrst density-modulus relationships for the ju- venile equine MC3 and P1. To evaluate whether indeed this relationship is di↵erent from those de- rived from older bone, we compared our juvenile 265 longitudinal MC3 data to that reported by Les et al [10]. Apparent CT density was converted to ash density ( ⇢ash(g/cm3), Eq. 1 [25]) and ash Figure 4: Density-modulus relationships for the (A) MC3 and (B) P1 samples in the longitudinal (orange) and transverse (blue) directions. Printed R2values are marginal R2values, which measures the variance in modulus that can be explained by CT density. In each of these relationships, subject (donor of the bone samples) was included as a random e↵ect. 5 density-modulus relationships were calculated. ⇢ash=(⇢CT+0.09)/1.14 (1) The root mean squared error (RMSE) was 347 270 MPa using Eadult,w h i c hi sa l m o s t8 0 %h i g h e r than the RMSE when using the model developed here for juvenile bone, Eyoung (191 MPa) (Fig. 5). Modulus values for adult MC3 trabecular bone in compression range from 2.09-3.65 GPa 275 [11, 12, 26], while median modulus for our juve- nile MC3 longitudinal samples was 609 ±297 MPa. It should be noted that over 80% of the sam- ples tested by Les et al. were cortical bone (al- though they reported no di↵erence in the density- 280 modulus relationship between cortical and trabec- ular bone [10]), while our samples were trabec- ular. Bone volume fraction in the distal MC3 condyles reportedly increases from approximately 32% in 1-2 month old horses to approximately 285 60% in horses greater than 6 years old [27]. As ap- parent CT density is a combined measure of bone volume fraction and tissue density, it is expected that apparent CT density changes with matura- tion. 290 While di↵erences in intrinsic properties be- tween immature and mature bone are the most likely explanations for why density-modulus rela- tionships di↵er with age, there are also method- ological di↵erences in the measurement of Young’s 295 modulus that are worth mentioning. We tested our samples using endcaps while Les et al. tested samples in compression using platens directly in contact with the bone sample, which has since been shown to result in an underestimation of 300 modulus by 20-40% [23] in trabecular bone sam- ples. 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 305 to be conﬁrmed. We found that density-modulus relationships in juvenile equine trabecular bone vary depending on anatomical location (MC3 vs P1) and load- ing direction (longitudinal vs transverse), which is 310 consistent with human data [18, 9]. The impact of not accounting for di↵erences between bones Figure 5: Density-modulus relationships for the MC3 comparing adult ( Eadult [10]) and juvenile ( Eyoung) samples. when predicting modulus from CT data is no- table. Applying the MC3 density-modulus rela- tionships to the P1 data results in an RMSPE 315 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%). 320 There are several reasons why the density- modulus relationships between bones may be dif- ferent. Biomechanically, the two bones are likely under di↵erent types of loads with the MC3 in combined bending and compression due to its long 325 slender structure and distribution of cortical ma- terial properties [28]. In contrast, the cuboidal P1 would be more resistant to bending. Varia- tion in surface strain modes between the distal MC3 (compression) and the P1 (shear) have been 330 reported [29, 30]. We also found that the density-modulus rela- tionships in the transverse direction had stronger predictions of modulus than the longitudinal di- rection (Fig. 4). The density-modulus relation- 335 ship of the P1 in the transverse direction ( EP1T) had the highest R2(0.86) and lowest percent error (RMSPE=23%) of all relationships investigated. 6 The density-modulus relationship of the MC3 in the transverse direction ( EMC3T)a l s oh a dah i g h 340 R2(0.77) but had the highest percent error (RM- SPE=45%) of all relationships, which was driven by increased variability of the modulus data be- tween a CT density of 0.2-0.5 g HA/cm3.F o r example, at a density of approximately 0.28 g 345 HA/cm3in the MC3 (Fig. 4A), the transverse modulus ranges from approximately 200-800 MPa (blue data). The sensitivity of the strength of modulus pre- dictions to orientation may be related to the na- 350 ture of the microstructure along each direction. Longitudinal cores tend to have more varied mi- crostructure along the length of the sample (Fig. 6A) when compared to the more compact mi- crostructure evident in transverse cores (Fig. 6B). 355 However, these microstructural di↵erences may be unique to juvenile animals, as Keyak et al. re- ported similar R2values in density-modulus rela- tionships for adult human proximal tibia bone in the longitudinal (0.84) and transverse directions 360 (anterior-posterior: 0.72; medial-lateral: 0.84) [9]. Augat et al. found approximately equivalent co- ecient of variation in modulus between the lon- gitudinal and both transverse directions of adult human trabecular bone in the spine, calcaneus, 365 proximal femur, and distal femur [31]. Further work is needed to conﬁrm whether mineralization rates are di↵erent between the MC3 and P1, as well as the inﬂuence of microstructure and tissue mineral density on the elastic modulus. 370 Although a driving motivator for this study was the ability to more accurately evaluate the me- chanical environment during exercise in juvenile horses, there are other applications. Finite ele- ment models based on CT data have been used 375 to assess fracture risk and represent an improved assessment of bone strength when compared to densitometric variables alone [32]. Finite element models can also be used as a pre-surgical evalua- tion tool to predict tissue response to certain ﬁx- 380 ation methods [33]. Use of the incorrect density- modulus relationship may result in incorrect as- sessments of bone strength or pre-surgical evalua- tion, therefore, our ﬁndings may have direct clin- ical implications. 385 Figure 6: MicroCT images of representative cores in the (A) longitudinal and (B) transverse directions of the MC3 and P1, from the same subject. All images are at nominal isotropic resolution of 144 µm. Scale bar is equivalent for all images. There are several limitations of this study. Transverse bone samples were only machined in the dorsal-palmar (“anterior-posterior”) direc- tion. The geometry of the MC3 and P1 led to challenges in excising cores in the medial-lateral 390 direction, and limited availability of intact juve- nile bones required we waste as little tissue as pos- sible. The implications of this may be mitigated by the fact that trabecular bone has been de- scribed as a transversely isotropic structure [34], 395 and modulus values in the anterior-posterior and medial-lateral directions often have a similar re- lationship with ash density [9] and bone volume fraction [35]. Despite the fact that we tested over 200 samples, we still encountered variability in 400 the density-modulus data that leaves over 30% of the modulus variability unexplained in the case of the longitudinal MC3 (Fig. 4A, EMC3L). Aside from the inﬂuence of microstructure on mechan- ical properties, some of this variability may be 405 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 sam- ples in the P1 (combining both orientations) were from subjects that were 4 weeks old. Sample size 410 in the current study does not allow us to develop statistically meaningful density-modulus relation- ships for each age of juvenile horses, and instead data were pooled to describe horses less than 1 7 year old. Nonetheless, all model predictions were 415 signiﬁcant and this work represents the ﬁrst large scale mechanical testing study in juvenile equine bone. Therefore, using rigorous imaging and experi- mental protocols, we have established orientation- 420 speciﬁc density-modulus relationships for the ju- venile 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 425 potential for bone adaptation. 5. Conﬂict of Interest None 6. Acknowledgements Samples were collected by AMM during 2018- 430 2022 from patients requiring euthanasia at the Large Animal Clinic of the Veterinary Teach- ing Hospital at UIUC. We are grateful for the technical support from Dr. Leilei Yin for his help with µCT scanning, and to the Visualiza- 435 tion Laboratory at the Beckman Institute for Ad- vanced Science and Technology for providing soft- ware used in this project. Funding was provided by USDA Hatch Funds (ILLU-888-964), Campus Research Board at UIUC, Grayson Jockey Club 440 Research Foundation, and a graduate fellowship from the Beckman Institute for Advanced Science and Technology to SGM. 7. Nomenclature CT - computed tomography 445 E-d e n s i t y - m o d u l u sr e l a t i o n s h i p MC3 - third metacarpal P1 - proximal phalanx Greek Letters µ-m i c r o 450 ⇢-d e n s i t y ,g / c m3Subscripts adult - data from Les et al. [10] ash - relating to amount of mineral, indicating ash density 455 CT - from clinical CT, indicating apparent CT density MC3-L - from the MC3 bone in the longitudinal direction MC3-T - from the MC3 bone in the transverse 460 direction P1-L - from the P1 bone in the longitudinal di- rection P1-T - from the P1 bone in the transverse direc- tion 465 young - juvenile MC3 data in the longitudinal di- rection