Finite element analysis (FEA) with models derived from computed tomography (CT) scans is potentially powerful as a translation al research tool because it can achieve what animal studies and cadaver biomechanics cannot – low-risk, non -invasive, objective assessment of outcomes in living humans who have actually experienced the injury or treatment being studied. The purpose of this study was to assess the va lidity of CT -based virtual mechanical testing with respect to physical biomechanical test s in a large animal model. Three different tibial osteotomy models were performed on 44 sheep . Data from 33 operated limbs and 20 intact limbs was retrospectively anal yzed. Radiographic union scoring was performed on the operated limbs and physical torsional test s were performed on all limbs . Morphometric measures and finite element (FE) models were developed from CT scans and virtual torsional tests were performed with four different material assignment techniques. In correlation analysis, m orphomet ric measures and radiographic scores were unreliable predictors of biomechanical rigidity, while the virtual torsion test results were strongly and significantly correlated with measured biomechanical test data, with high absolute agreement . Overall, the results validated the use of virtual mechanical testing as a reliable in vivo assess ment of structural bone healing. T his method is readily translatable to clinical evaluation for noninvasive assessment of the healing progress of fractures with minimal risk. Clinical significance: Virtual mechanical testing can be used to reliably and non -invasively assess the rigidity of a healing fracture using clinical -resolution CT scans and that this measure is superior to morphometric and radiographic measures . Keywords: Finite element analysis; Computed tomography, Tibial shaft f racture; Ovine osteotomy; 3 Introduction Finite element analysis (FEA) with models derived from computed tomography (CT) scans is increasing in prevalence as a research and clinical diagnostic tool in orthopaedics .1–8 This technique is able to accurately capture both the complex geometries and localized material acquired by the CT scan to build representative animal -specific or patient -specific models .9,10 FEA can then be used to simulate different loading scenarios and assess the structural biomechanics of the anatomy of interest . With animal studies and cadaver biomechanics, structural outcomes can be directly measu red. However, outcome assessments such as benchtop biomechanical testing and histology cannot be performed in clinical settings with human patients. To that end, image -based models are potentially powerful as a clinically translatable research tool because they can achieve what animal studies and cadaver biomechanics cannot – low-risk, non -invasive , objective assessment of outcomes in living humans who have actually experienced the injury or treatment being studied . For example, we rece ntly developed CT -derived computational methods to calculate the virtual torsional rigidity (VTR) of human tibial fractures after 12 weeks of healing .9 In a pilot study, we showed that this measure is strongly correlated with clinical healing outcomes such as time to clinical union.10 Furthermore, integrating image -based models in preclinical in vivo studies may reduce the required number of experimental animals, in accordance with 3R principles, by allowing non -invasive follow -up assessments of structural healing kinetics. One barrier to widespread clinical use of patient -specific modeling techniques is that unlike in preclinical animal models, the quantitative predictions from these models cannot be easily validated due to a lack of physical biomechanical testing data in hu mans . Previous investigators have tried to bridge the gap between preclinical and clinical research techniques by using 4 biomechanical testing of small animal osteotomy models to validate clinical assessments such as the radiographic union score (RUS or RUS T when applied to the tibia). Recently, two independently conducted murine studies compared radiographic assessments with biomechanical data and CT -derived callus morphometry measures .11,12 When interpreted together, these studies suggest that radiographic sc oring is less reliable (lower inter -observer consistency) and less predictive of biomechanics (weaker correlation) at early timepoints in normal healing condition and with compromised fracture healing conditions compared to later timepoints in normal healing conditions. One previous study examined the repeatability of radiographic scoring when applied in well -healed ovine tibial osteotomies, but did not report correlations with the results of physical biomechanica l testing.13 These gaps in our ability to assess fracture union quantitatively across the spectrum of healing could be resolved by adoption of CT -based virtual biomechanics, but the ne cessary validation data showing that the virtual tests replicate the physical results have not yet been presented. Accordingly, the purpose of this investigation was to assess the validity of CT -based virtual mechanical testing with respect to physical bio mechanical testing in a large animal model. For this purpose, we used a series of standardized ovine (sheep) tibia osteotomies stabilized with internal fixation for which both CT scans and postmortem biomechanical testing data were available and which repr esented a wide spectrum of healing responses. The hypothesis of this study was that the results of physical and virtual biomechanical testing would be strongly correlated. Methods Specimen Information: 5 Forty -four adult female Swiss alpine sheep (2-3 years old, weighing 59 – 87 kg) were part of two previously completed research studies with three different tibial osteotomy models (Figure 1) stabilized by medial plating . In total, there were 33 operated limbs and 20 intact control limbs that had both mechanical test data and CT scans available for analysis in this study . Taken together, these animals comprised three experimental datasets across a wide spectrum of healing responses . Dataset 1 consisted of data from seven animals with a 3 mm gap defect stabilized with a 12-hole stainless steel plate ( broad straight veterinary 3.5 mm locking compression plate (LCP) , 159 mm in length, with 3.5 mm bicortical screws; DePuy Synthes® ). Dataset 2 consisted of data from 18 animals with a 3 mm gap defect stabilized with a six-hole titanium plate (broad 4.5/5.0 mm LCP, 115.8 mm in length, with 5 mm bicortical screws; DePuy Synthes® ). Dataset 3 consisted of data from eight animals with a 17 mm defect augmented with autografts and stabilized with a 13 -hole stainless steel plate (broad straight veterinary 3.5 mm LCP, 172 mm in length, with 3.5 mm bicortical screws; DePuy Synthes® ). The 3 mm defect models (Datasets 1 and 2) represented a non -critical size defect capable of spontaneously healing . Sheep in Datasets 1 and 2 were sacrificed at 9 weeks. The 17 mm graft model (Dataset 3) represented a critical size defect that would not spontaneously heal without autograft augmentation . Sheep in Dataset 3 were sacrificed 12 weeks after surgery. All animal experiments were conducted according to the Swiss laws of animal protection and welfare and authorized by the local governmental veterinary authorities (license numbers ZH071/17 and ZH183/17) . Radiographic Union Scoring : Plain film radiographs w ere taken postmortem at the time of sacrifice in three different planes: anteroposterior (AP), mediolateral (ML), and an angled plane between the AP and ML viewed from the cranial aspect of the limb. Radiographs were evaluated by two independent, 6 board cer tified expert reviewers (Prof. Dr. med. vet. Brigitte von Rechenberg , Dipl. ECVS and Prof. Dr. med. vet. Mark Flückiger , Dipl. ECVDI ). Semiquantitative radiographic union scores were assigned following a scoring system based in part on the modified radiographic union score for tibial fractures (mRUST ) method.14–16 Our scoring approach included the callus bridging assessment of clinical mRUST , except that due to medial plate fixation , only the lateral callus was scored from the AP view due to the presence of the medial plate . Half-point scores were allowed for greater granularity on the uniform osteotomies . The resulting unicort ical score is hereafter referred to as lateral callus granular mRUST [score range 1 - 4 in 0.5 increments ]. The callus bridging scores were performed on the AP radiograph . Additionally, the following scoring criteria we included for callus maturity, which are not typically assessed in clinical mRUST: cortical callus formation (AP projection ) [0-4], cis-cortex callus formation (AP projection ) [0-4], trans -cortex callus formation (AP projection ) [0-4], cranial cortical gap (ML projection ) [0-4], caudal cortical gap (angled projection ) [0-4], and callus opacity scores (AP projection ) [0-3]. All scoring components were summed, resulting in a comprehensive radiographic union score (cRUS) on an interval range [1-27] with higher scores representing more advanced healing . A scoring breakdown can be seen in the supplementary digital content. Micro -Computed Tomography (μCT) Scanning: After animal sacrifice, tibiae were excised, stripped of soft tissue , wrapped in saline -soaked gauze, and µCT scanned using an XtremeCT II Micro -CT scanner (Scanco Medical AG, Bruettisellen, Switzerland) with an X -ray voltage of 68 kVp and X-ray current of 1 470 µA. The resulting scans had an isotropic resolution of 60.7 µm . A phantom (Scanco KP70 phantom, QRM, Moehrendorf, Germany) was scan ned in the same scanner at identical settings allowing 7 for conversion from native Hounsfield Units ( HU) to calibrated mineral density ( ρQCT, mgHA/cm3). Biomechanical Data Physical torsion tests were performed on all included samples using a custom -made fixture on an Instron E10000 electro dynamic testing machine ( Instron, Massachusetts, US ). Axial loads and torques were measured with a calibrated load cell (±10 kN / ±100 N -m). Each tibia was prepared by stripping the periosteum on both ends and fixing the tibia in the test frame using Beracryl embedding medium . Periosteum was not stripped in the region of the callus, only at the embedded ends. To minimize bone movement in th e Beracryl, four adjunctive screws were inserted perpendicular to the mechanical axis of the bone in the distal epiphyses. Potting depth was adjusted to ensure a minimal distance of 10 mm between the Beracryl and the nearest screw hole from the fracture pl ate. The spacing between the proximal and distal Beracryl surfaces after embedding was recorded and ranged between 140 mm and 160 mm. While potting, the diaphysis was kept moist using saline -soaked gauze. Mechanical tests were performed by quasi -statically preloading the limb with 5 N axially, which was held for the entire test, and then applying internal rotation at 5 °/min. Biomechanical torsional rigidity was calculated as a linear regression of the loading curve between 6 and 10 Nm multiplied by the gauge length, or the distance between the surfaces of the two Beracryl pots, of the test specimen during testing. Scan processing and 3D model construction: CT scans were processed using the Mimics Innovation Suite (Materialise, Leuven, Belgium) following a similar work flow as Schwarzenberg et al.9 To better replicate the c linical scans used in the previous study and ensure future translatability , all of the µCT scans were 8 initially down -sampled to an isotropic resolution of 400 µm prior to segmentation . Density threshold values of 400 -2500 HU and 2500 -4000 HU were chosen to initially segment the callus and cortical bone respectively (Figure 2c-d). Proximal and distal cut -planes were selected just interior to the most proximal and distal screw holes, resulting in two planar cortical surfaces for FE boundary condition (BC) application. Following preliminary threshold -based segmentation of callus and cortical bone, all models were screened for formation of high -density tissue at the cortical -callus boundary, which can occur in samples with more advanced healing. To accomplish this, each callus mask was reviewed in a slice -by-slice manner throughout the CT -stack. In each slice, voxels of high - density new tissue that fell outside the contours of the original callus mask were assigned to a new tissue mask region. This step ensured that the callus morphometric measures were accurate and that the cortical bone fragments represented exclusively the original bone contour and none of the new tissue growth. After the mask for the new tissue was segmented , it was subtracted from the density -based cortical bone mask to create the true pre-existing cortical bone mask (Figure 2 e-f). These regions are referred to as callus and bone within this work . Final geometries were visually inspected by a surgeo n to ensure the virtual models were anatomically representative of the biological samples . The final models for all operated limbs can be seen in Figure 3. After the callus and bone masks were accurately segmented , the masks were united into a single surface model for preparation of volumetric discretization. The united models were wrapped with a gap closing distance of 1 mm and a smallest detail of 0.5 mm to produce cohesive surfaces . Finite element model creation : 9 The united surface model s were exported to a dedicated toolkit in the Mimics Innovation Suite, 3 -Matic, where the surface models were smoothed to further reduce CT scan noise and to fill any small defects. Next a triangular mesh was appli ed to the surface of each model with maximum edge lengths of 0.4 mm and a linear tetrahedral volumetric mesh was applied to the body of each model with maximum interior edge lengths of 0.875 mm. These parameters are based on a previously published mesh conve rgence study conducted with human fracture callus models9. Elementwise material properties were initially applied to the FE mesh es (Figure 2 h) by three non-species -specific methods : one linear and two power law material assignments taken from the literature (Figure 4 ). First, a direct conversion from HU to elastic modulus was achieved using a linear scaling law .17 For the two power law material assignment s, the native HU density measure of the CT scan was first converted to calibrated mineral density (ρ QCT) using the phantom and then to ash and apparent density using equations from Schileo et al.18 Once voxel densit y values were in ash and apparent densities, material scaling law equations were used from Keyak et al.19 and Morgan et al.20 respectively to calculate elastic modulus. These two power law equations were chosen to represent the lower and higher end of elastic modulus assignment from two literature review s.21,22 A species -specific scaling law for ovine tibial cortical bone was not identified f rom a literature search and was subsequently developed in this study. Virtual torsion testing In preclinical studies, torsion testing is the preferred mode of assessment due to its robustness to variation in setup orientation . Structural FE simulations rep licating this torsional test were carried out in ANSYS 17.2 (Canonsburg, Pennsylvania). The boundary conditions were rigid fixation on the distal end and 1 degree of applied rotation on the proximal end, with twisting 10 about the aligned mechanical axis of t he tibia. Additionally, using the same multi -point constraint (MPC) contact of the applied axial rotation, the proximal end was fixed in the anterior - posterior and medial -lateral translational direction s to best replicate the physical test setup (Figure 2 i). For each model, the virtual torsional rigidity (VTR) was calculated as the moment reaction from the applied loading ( 𝑀) multiplied by the working length of the test segment ( 𝐿) divided by the applied angle of twist ( 𝜙): 𝑉𝑇𝑅=𝑀𝐿/𝜙 [Nm2/°]. Material Optimization An optimization technique similar to Eberle et al.23 was performed on the N = 20 intact tibiae to find a material assignment law of the form 𝐸=𝑎𝜌𝑄𝐶𝑇𝑏 that best fit the data. In brief, value ranges for a (5,000 ≤ a ≤ 20,000) and b (1 ≤ b ≤ 3) were chosen from the same literature review as above21,22 to represent the current published data across multiple anatomic sites . Sixteen design points were chosen to span the parameter space and the torsional rigidity of each intact model was calculated by assigning all elementwise material properties within all N = 20 models using the coefficients a and b at those points. A root mean square error (RMSE) was calculated comparing each design point (coefficients a and b) to the biomechanical data and a 2nd order response surface was fit to the a, b, and RMSE values. The minimum predicted RMSE of the response surface was calculated by 120,801 sample points in the design space , defined by a resolution of 25 (𝑀𝑃𝑎𝑐𝑚3 𝑔𝐻𝐴) in the a dimension and 0.01 (-) in the b dimension. After the initial response surface minimum was calculated, three additional refinement design points were iteratively computed from the surface minimum. Finally, the minimum a and b values for the response surface with 19 design points (16 preliminary and 3 refinement) was used to assign elementwise material properties in all N = 33 operated limb models and N = 20 intact limb 11 models . The resulting virtual torsional rigidities ( VTRs ) were calculated for further statistical analysis . Statistical Analysis Descriptive statistics were generated using Microsoft Excel and MATLAB (R201 6a- R2019 a, The MathWorks, Inc., Natick, Massachusetts). Additional statistical analys es were conducted in SPSS 25 (IBM Corp., Armonk, NY). Pearson correlations were performed to assess the strength of the linear association between various outcome measures and the measured biomechanical rigidity from physical testing . Correlation coefficients were interpreted as follows24,25: weak R2 ≤ 0.3, moderate 0.3 < R2 < 0.6, or strong R2 ≥ 0.6. One-way repeated - measure s ANOVAs were computed to determine group differences between the biomechanical rigidity and the simulated VTR data. Additionally, root mean squared error ( RMSE ) analys es were performed to evaluate the relative performance of the four material property assignment technique s (three non -species -specific, one ovine -optimized) . All reported values are medians and interquartile ranges (IQR) unless otherwise reported and statistical significance was at p < 0.05. Results Structural Data: Biomechanical Testing The results of physical biomechanical testing for all samples are shown in Figure 5 and Table 1. Dataset 3 (17 mm defect with graft) had the lowest torsional rigidity, as expected in this delayed healing model. Between -groups comparisons were not performed between Datasets 1, 2 and 3 because these experiments were not conducted concurrently with intent to compare outcomes. Instead, t hese torsional rigidity results were compiled into a single osteotomy dataset , 12 representing a range of healing responses, and were subsequently used as a benchmark to assessing the reliability of various outcomes measures (morphometric, radiographic, and virtual biomechanics) with respect to the measured physical properties of the specimens. Morphometric Data Summary morphometric data for the three osteotomy datasets individually and for all samples together is presented in Table 2. Callus volume was weakly -to-moderately and significantly correlated with biomechanical torsional rigidity in Datasets 2 and 3 and for the 33 compiled osteotomies overall. Callus density was weakly and non -significantly associated with biomechanical torsional rigidity in all datasets and for the 33 compiled osteotomies overall. Radiographic Data Radiographic scoring data can be seen in Table 3. The lateral callus granular mRUST score was moderately and significantly correlated with biomechanical torsional rigidity in Dataset 1 only, but weak ly and non -significant ly correlated with the other two datasets and the 33 compiled osteotomies overall. Similarly, the comprehensive radiographic union score (cRUS) was strongly and significantly correlated with biomechanical torsional rigidity in Dataset 1 only, with weak and non -significant correlations for the oth er two datasets and the 33 compiled osteotomies overall. Structural Data: V irtual Torsional Rigidity (V TR) Virtual torsional rigidity (VTR) data from all simulations is shown in Figure 5 and Table 4. As described above, each of the N = 33 osteotomy limbs and N = 20 intact limbs were subjected to virtual testing with four different material assignment approaches: three non - species -specific scaling laws (Synder linear, Morgan power law, and Keyak power law) and one 13 optimized sp ecies -specific ovine power law. Across all osteotomy datasets , the compiled osteotomies, and the intact limbs, regardless of material assignment law, the virtual and physical biomechanical tests were moderately -to-strongly and significantly correlated. Uns urprisingly, the VTRs resulting from the optimized ovine -specific material law performed as well or better than the non -optimized non -species -specific material assignment approaches. Figure 5 compares the predicted VTR values for the compiled osteotomies and the intact limbs to the physical biomechanical torsional rigidity. Significant differences are indicated with respect to physical biomechanical torsional rigidity. The non -species -specific material assignment laws significantly over -predicted rigidity , despite their high correlation with the biomechanical test results. T he RMSE analysis showed that in both the fractured and intact models, the Keyak power law material assignment method had the lowest relative error of the non-species -specific material assignment laws (Table 5). Discussion The goal in many preclinical orthopaedic studies is to measure the healing progress of the fractured bone with a surrogate measure such as torsional stiffness or rigidity. To accomplish this, animals need to be sacrificed and bones excised for physical testing, which cannot be done in clinical studies or at interim timepoints in an animal model . Currently, morphometric measure ment of callus formation (volume and density) is a commonly used technique that can b e performed in vivo to evaluate fracture healing progress prior to animal sacrifice .26,27 Furthermore, the clinical gold -standard radiographic scoring systems such as RUST and mRUST a re used in both preclinical and clinical settings to track healing progress . In this study , our method of virtual torsion testing outperformed all of these measures when predicting the biomechanical rigidity of the b one. Virtual torsion testing with the species -specific material model also 14 produced strong and statistically significant correlations and low error between the measured and simulated biomechanical properties across a wide range of healing responses . The study results also indicate that morphometric measures of callus formation should be interpreted with caution. Callus volume and density are often considered surrogate measures of healing progress o r even indirect indicators of biomechanical properties . In this context, c allus volume is analogous to the amount of new tissue visible on plain fil m radiographs , while callus density can help infer tissue maturity and structural integrity . Both are interpreted as signs of healing .28 In this study, callus volume had a moderate co rrelation with measured torsional rigidity in Dataset 1 and Dataset 3 and weak correlations with torsional rigidity in Dataset 2. Dataset 3 was designed to represent a critical size defect which is more indicative of a delayed healing case which may have more dependency on callus volume than a non -critical size defect model. Furthermore, this dependency on callus volume may be attributed to Dataset 3 being a delayed healing model in which new bon e is still expected to significantly remodel and condense after the analyzed endpoint of 12 weeks . However, within all datasets , callus density showed little association with biomechanics , which could be because all callus present was of similar maturity . In a preclinical model with more widely varying timepoints, callus density could be a more important indicator of healing progress . In our controlled osteotomy models the coefficients of variation for callus volume (36%, 38% and 24%) were much higher than for callus density (10%, 8% and 7%). This is consistent with callus volume and density patterns observed in humans at 12 weeks10. Notably, the plates used for fixation were not the same across datasets , which could explain some of the differences in healing response, although as previously stated, between -groups comparison of the fixators was not the objective of this investigation . In fact, by design we aimed to show the validity of our hypothesis across various fixation 15 techniques and defect sizes. The observed low correlation of morphometric measures with outcome parameters seen in this study are consistent with findings from our clinical pilot studies , which used the same virtual mechanical testing technique and showed that morphometric parameters are not strongly correlated with simulated VTR or with clinical outcomes such as time to union .10 Additionally , the radiograp hic scoring results do not provide a consistent picture. T he radiographic scores in Dataset 1 had much stronger correlations with biomechanical rigidity than the other two models . The callus in Dataset 1 also tended to be the largest (Table 2). This could suggest that visual radiographic scoring is more reliable when there is more callus present. If true, this would be consistent with murine data showing that radiograph ic scoring is less reliable and less predictive of biomechanics at early timepoints and with compromised fracture healing compared with later timepoints in normal healing conditions29. Overall, t his study validated the robustness of virtual mechanical testing as a n in vivo method for assessing structural bone healing . With the optimized species -specific material assignment law for ovine tibial cortical bone, w e found a strong and statistically significant correlation between virtual and physical torsional rigidity across all sample subgroups , which was superior to both radiographic and mo rphometric assessment s from CT . Even with non -species -specific material models, the virtual mechanical test results were still moderately -to-strongly and significantly correlated with physical biomechanical testing. This finding is particularly notable because it indicates that virtual mechanical testing may be a powerful in vivo assessment tool, even when a robust species - and site -specific material assignment law is not availabl e. This study had a few additional minor limitations that may be restricting the predictive power and correlations between measures. First, as in all physical biomechanical torsion test s, there 16 were small potting artifacts where the bone was able to sligh tly deform within the PMMA pot. This could reduce the measure d stiffness and rigidity of the biomechanical data, raising the relative error between physical and simulated data. Additionally, while higher resolution µCT scans were acquired , the scans were d own sampled to replicate the resolution of clinical scanners . Starting with a less resolute data set eliminates the possibility of more sophisticated structural material modeling, such as the use of fabric tensors, and could have also introduce d mask boundary definition errors in the virtual models. However, µCT scan data is not available in a clinical setting or at interim in vivo timepoints in large animal models, so the clinical resolution chosen here has clear translational relevance. Finally, it is important to note that the optimized scaling law we developed for the ovine tibial diaphysis using the intact limbs may not accurately model the properties of callus tissue at all stages of maturity. At the timepoint our scans were taken , fracture c allus consist ed of fibrous tissue, some cartilage, and woven bone, which is not as organized or mineralized as cortical bone. During coupled remodeling, callus microstructure, composition, and mechanics change adaptively30, and the modeling approach we are using my not be adequate to capture all of those processes. Optimization of new material property scaling laws for combined cortical -callus structures should be a target for continuing research and would undoubtedly reduce the RMSE of the virtual prediction. Conclusions Image -based structural mechanics modeling from CT scans is a n objective, quantitative method for assessing fracture healing that can closely replicate physical biomechanical testing . Using the methods described herein , we validated a technique for virtual torsional testing using ovine tibial osteotomy fracture models representing a wide range of healing scenarios . The 17 virtual mechanical testing results were better predictors of trends in measured biomechanical properties than radiographic scoring or morphometric analysis of callus from CT scans in these models. This method is readily translatable to preclinical or clinical settings for noninvasive assessment of the structural healing progress of fractured long bones with minimal risk. Acknowledgements The preclinical stud ies from which this data was obtained we re funded by the Johnson & Johnson Family of Companies and Biomech Innovations . The authors wish to thank Beat Lechmann (Johnson & Johnson Family of Companies) and Stefano Brianza (Biomech Innovations AG) for agreeing to grant access to the ovine stud ies data for our analyses. The authors also wish to thank Beat Lechmann (Johnson & Johnson Family of Companies) for Figure 1 graphics and Prof. Mark Flückiger for performing radiographic assessments together with BvR . The authors want to thank the National Sc ience Foundation for supporting Peter Schwarzenberg with a n IIE-GIRE travel grant. BvR is a paid consultant and receives research support as a principal investigator from Johnson & Johnson Family of Companies. HLD has stock or stock options, receives resea rch or institutional support as a principal investigator, and receives financial or material support from OrthoXel. Additionally, HLD serves on a committee in the Orthopaedic Research Society (ORS).