Bone fractures commonly repair by form ing a bridging structure called callus, which begins as soft tissue and gradually ossifies to restore the rigidity of the bone. In large animals and humans, this process typically requires only a few months to complete . Occasionally bone s can fail to heal, resulting in nonunion s that are difficult to diagnose and treat, partly because early detection of poor healing remains challenging. We present a translational technique fo r measuring what really matters in fracture healing – the structural integrity of the bone – using image -based virtual mechanical testing. This work show s that physical mechanical tests of osteotomized ovine tibiae can be reliably recapitulated in a simula tion environment only when the soft -hard mechanical duality of the healing zone is accurately captured by the constitutive material model. A large -scale optimization analysis was performed on thousands of 3D finite element models derived from computed tomo graphy (CT) scans of 33 osteotomized sheep. This identified a piecewise material model that successfully replicated the postmortem torsion testing data by differentiating between soft and hard callus based on material density alone. The results suggest tha t most of the structural integrity of a healing bone is conferred by an internal architecture of mineralized hard callus that is supported by interstitial soft tissue. We conclude that with appropriate material modeling, virtual mechanical testing is a rel iable surrogate for physical testing and that this technique has high translational potential for diagnostic testing of nonunion. 3 Significance Statement Broken bones often heal by forming new soft tissue around the injury site that gradually transit ions into hard bone. This work show s how this dual nature of the healing zone , as both a soft and hard material , determines the mechanical rigidity of the whole bone. Medical imaging data was used to build 3D models of healing sheep bones and analyze thous ands of variations of the soft -hard material behavior . Ultimately, w e identif ied a method that reliably replicate d physical mechanical tests of the bones in a virtual simulation environment. Our findings show the potential for using imag e-based mechanical simulations to accurately measure bone healing and potentially improve clinical care of patients who have poorly healing fractures. 4 Introduction Long -bone fractures are among the most common and painful trauma -associated injuries of the musculoskelet al system (1, 2) . Under conditions of relative stability, secondary fracture healing progresses by both intramembranous (periosteal) and endochondral bone formation to produce a structurally heterogeneous, functionally graded external bridging callus (3). Most fractures heal within a few months, but in some patients, the process fails, resulting in a nonunion (4). These patients have difficulty returning to work, report worse health -related quality of life than patients with AIDS, stroke, or type -I diabetes , and would give up 40% of their remaining lifespan to regain their health (5–7). Nonunion s remain a pervasive clinical problem in part because diagnosis currently requires subjective assessment of X -rays over long periods of watchful waiting, typically at least 6 -9 months (8). Unlike in anim al models, clinical interpretation of the radiographic progress of fracture healing relies on a surgeon’s judgement of what can be seen, rather than an objective measurement of the structural mechanics of the bone (9). To address this clinical need, we previously develop ed a non -invasive, imaging -based virtual biomechanical test for assessing the progress of fracture healing. First, w e developed a methodology for processing computed tomography (CT) scans to develop subject -specific 3D finite element models of healing bones (10, 11) . In a recent ovine study , we demonstrated that a virtual torsion test outperforms subjective methods like radiographic scoring for predicting the progress of healing and t hat it is a reliable surrogate for postmortem physical torsion testing in intact tibiae (12). In clinical pilot testing, out virtual torsion tests have been able to detect clinically significant delayed fracture healing in patients with comorbidities such as smoking and diabetes (13, 14) . 5 The key to success in virtual mechanical testing of bone mechanics is the use of density - dependent scaling laws for converting the radiodensity (gray value) in each voxel of the scan to a Young’s modulus , 𝐸, (15–17): 𝐸=𝑎𝜌𝑏 [1] where the coefficients 𝑎 and 𝑏 typically describe a linear or weak power -law relation and the radiodensity , 𝜌, is often expressed in either Hounsfield units [HU] or a calibrated bone mineral density using a radiological phantom . It is well understood that bone density is positive ly correlated with mechanical properties including Young’s modulus (18). In subject -specific fini te element models, elementwise material properties are mapped from the underlying scan data, such that regions of bone with lower density are assigned lower modulus values . This technique is highly effective for mimicking physical testing of intact bones b ut has a clear limitation when extended to be used in specimens with callus. Fracture callus is highly heterogeneous, as depicted in Fig. 1. This naturally adaptive structure develops in response to biological and mechanical cues and is comprised of hard c allus (woven bone, newly mineralized or approaching the density of the intact cortex) and soft callus (fibrous and/or cartilaginous interstitial tissues) (19, 20) . Our previous ovine study found that simply exten ding a density -modulus scaling law developed for cortical bone to include regions of callus produced strong correlation s between virtual and physical biomechanical tests, but that the virtual tests over -predicted the measured torsional rigidity in osteotom ized specimens by an average of 58% (12). This over -prediction of rigidity when callus is present suggests that a constitutive material model derived for cortical bone does not completely capture the mechanical behavior of callus and that more work is needed to virtually replicate the in vivo biomechanics of healing fractures. 6 Accordingly, t he purpose of this investigation was to develop a new constitutive mechanics model for fracture callus that captures its dual nature —hard and soft —and to test this method for assigning material properties in finite ele ment models of healing long bones . The hypothesis was that the virtual tests performed using t he new dual-zone material assignment law for callus would outperform the virtual tests that use current single -equation methods when validated to physical biomech anical data. Results Effect of Dual -Zone Material Model on Torsional Rigidity Selected results from the material model optimization procedure are shown in Fig. 2 . The reference single -zone material model (Equation 3) did not differentiate between soft and hard callus within the healing zone and it overpredict ed the virtual torsional rigidity ( VTR ) of the whole structure (Fig. 2A , up arrows) . In the dual -zone material model (Equation 4 ), any element with a calibrated radiodensity ( 𝜌𝑄𝐶𝑇) below the so ft callus cutoff , 𝜌𝑐𝑢𝑡, was assigned the soft callus modulus value ( 𝐸𝑠𝑐 = 0.5, 5, 50, or 500 MPa). Within the optimization space, one limiting case was for 𝜌𝑐𝑢𝑡 = 0 mgHA/cm3 for all 𝐸𝑠𝑐 levels. In this case , the dual -zone material model performance was identical to the single -zone material model, with no elements assigned the soft callus modulus value and all elements treated as bone of varying density . Increasing 𝜌𝑐𝑢𝑡 progressively increased the number of elements within th e healing zone that were assigned the soft callus modulus. As 𝜌𝑐𝑢𝑡 was increased, the virtual models transitioned from over -predicting to under -predicting rigidity compared to the biomechanical testing data. In the extreme case, increasing 𝜌𝑐𝑢𝑡 into the range for cortical bone assigned the soft callus modulus to all elements in the model and resulted in very low torsional rigidity for the structure ( Fig. 2A , down arrows). 7 The optimization process was successful in identifying a value for 𝜌𝑐𝑢𝑡 at each level of 𝐸𝑠𝑐 that minimized the difference in torsional rigidity between the virtual and physical biomechanical tests ( Fig. 2B -D). At the optimized combination of soft callus cutoff and modulus, the predicted VTR converged with biomech anical torsional rigidity (GJ) and the RMSE dropped to its minimum value ( Fig. 2B/C). An interactive web application that illustrates the effect of changing 𝜌𝑐𝑢𝑡 on the performance of the dual -zone material models can be found in SI Appendix , Fig. S 1. Optimization Results Summary results of the optimization can be found in Fig. 3 and Table 1. The optimized density cutoffs for each soft callus modulus value (0.5, 5, 50, and 500 MPa) were determined to be 610, 625, 665, and 825 mgHA/cm3 respectively . The performance of the optimized dual -zone models at all four of the candidate soft callus modulus values is visualized in Fig. 3 . As previously reported, the virtual mechanical tests with single -zone material assignment reliably recapitulated the physical mechanical tests for the intact ovine tibiae, with a strong correlation and good absolute agreement ( Fig. 3A; R2 = 0.697 , RMSE = 0.10 6) (12). Applying the single - zone material model in the operated limbs also produced a strong correlation, but with VTR clearly over -predicting GJ (Fig. 3B; R2 = 0.638 , RMSE = 0.451). The absolute agreement did not meet our defined criteria for success (RMSE < datum standard deviation (SD); see Table 2 ). In contrast, at all tested levels of soft callus modulus, the optimized dual -zone material models produced moderate -to-strong correlations between VTR and GJ with absolute agreement that met our defined success criteria for material model validation (all R2 ≥ 0.558 and RMSE ≤ 0.248 versus datum SD = 0.268 ; Fig. 3C/D and Table s 1 and 2 ). In both the Training and Testing groups , the optimized dual -zone model with 𝐸𝑠𝑐 = 500 MPa had the lowest RMSE at 0.14 8 and 0.176, respectively. The model with 𝐸𝑠𝑐 = 50 MPa had the second -lowest Training and Testing 8 RMSE s at 0.160 and 0.203, respectively . Section contour plots showing a comparison of the material assignment distribution within the callus for the single -zone and optimized dual -zone material models ha ve been included for representative animal s in SI Appendix , Fig. S2. Across both the Training and Testing groups , VTR predictions with the optimized dual - zone materials showed better agreement with GJ compared to the single -zone material model (Fig. 3E/F). All four of the candidate dual-zone material models had interquartile ranges that overlapped the interquartile range of the biomechanical data, while the single -zone law clearly overpredict ed the stiffness. ANOVA post-hoc testing showed statistically sig nificant differences between the single -zone VTR and biomechanical GJ for both the Training and Testing groups (both p < 0.001, Fig. 3E/F ). There were no significant differences between biomechanical GJ and VTR for any dual-zone material model in the Training group . In the Testing group , there were no significant differences between GJ and VTR in the dual -zone models for soft callus modulus values of 5 MPa and higher; the only significant difference relative to the biomechanical test group was at the lowest soft callu s modulus level ( 0.5 MPa ; p = 0.007 , Fig. 3F). Physiological Performance Evaluation: Material Properties, Strain , and Histology For the Training group , all four of the optimized dual -zone material models produced correlation coefficients that were higher than tha t achieved with the single -zone material model (𝑅2 ≥ 0.724 vs. 𝑅2 = 0.638), with better absolute agreement with the physical testing (RMSE ≤ 0.185 vs. RMSE = 0.451 ). Correlation coefficients were slightly lower in the Testing group , but absolute ag reement between VTR and GJ was still clearly improved relative to the single -zone material model. Except for the lowest soft callus modulus level ( 𝐸𝑠𝑐 = 0.5 MPa), which was significantly different from the physical testing in the Testing group , the r esults did not provide a 9 statistical justification for making a recommendation on the best -performing dual -zone material model ( Fig. 3E/F and Table 1 ). To gain additional insights on possible differences in the structural model performance between these groups , contour plots of equivalent (von Mises) strain and Young’s modulus in the callus region were qualitatively compared to the histological sections. Example contour plots and section images are shown in Fig. 4 for one representative animal and each of the four candidate dual-zone material models (optimized combinations of 𝐸𝑠𝑐 and 𝜌𝑐𝑢𝑡). For low values of the density cutoff , 𝜌𝑐𝑢𝑡, the callus region contained a large proportion of unmodified elements that were assigned material propert ies using the cortical bone scaling law (Equation 3 ), resulting in a stiffer structure with low strains throughout . Strain contour plots showed that increasing 𝜌𝑐𝑢𝑡 was also associated with increasing concentrations of strain in and around the fract ure gap. At high 𝜌𝑐𝑢𝑡, the strain field became non -physical with poor correspondence to histological patterns . At the optimized combination s of 𝐸𝑠𝑐 and 𝜌𝑐𝑢𝑡, the assigned mechanical properties and strain contours were qualitatively well -matched with the localization of less -mature callus in the histological sections. Additional side -by-side comparisons of strain and histological slices can be found in SI Appendix , Fig. S3. Discussion The d evelopment of a clinically relevant, reliable, scalable, non-invasive test to assess bone healing could have profound impacts on the care of nonunion fractures. Detecting problems in bone healing earlier would enable earlier intervention when warranted, and potentially better outcomes for the most at -risk patients. Image -based virtual mechanical testing shows promise for addressing this need, but one key barrier to clinical translation of this technique is the problem of 10 material modeling in fracture callus. Callus is an inherently transient material t hat persists in a living organism for only a few months before being remodeled. Unlike native bone, the mechanical properties of callus cannot be studied in human cadavers, so our understanding of the material mechanics and structural organization of this tissue must be strongly informed by observations from large animals. The results of this study clearly demonstrate that t he material properties of the tissue s in the callus region for ovine long bones are not monolithic and should therefore not be treated as such. This idea has strong resonance with h istomorphometric data showing subzones of different composition within fracture callu s and gradients in mineralization that develop in both space and time (21–23). Vetter et al. (2010) performed quantitative tissue type mapping on sections from osteotomized ovine tibiae across the stages of early healing and found dramatic heterogeneity of density in regions where bone was form ing, with intermixing of bony and fibrous tissues and a moving front of bone formation (24). This agrees with nanoindentation mapping data that revealed local structural heterogeneity within small regions of callus and evi dence that calcium deposition correlates with local tissue modulus (25). The dual-zone or piecewise material model adopted in this study helps to capture callus heterogeneity in two ways . First, it recognizes that where the callus contains bone, this tiss ue should be locally modeled as bone. To achieve this, we extend ed our previously validated scaling law for ovine cortical bone to include the radiodensity range of woven bone within the callus . Second, the dual zone material model acknowledges that in ear ly or slow healing, a substantial portion of the healing zone may be comprised of fibrous and cartilaginous soft tissues. The results of this study confirm that from a structural mechanics perspective, soft tissues within the callus are not functioning sim ply as low - density bone . This can be inferred because using a bone -derived radiodensity scaling law to 11 assign their mechanical properties leads to a systematic over -prediction of limb rigidity. In contrast, by dropping the assigned soft tissue modulus for low-density regions within the callus, strong correlations between the virtual and physical mechanical tests were achieved with greatly improved absolute agreement compared to a single -zone material model approach. Importantly, these results were achieved using a purely objective, automatic density - based thresholding technique for identifying soft callus that d oes not require histological analysis, contrast enhancement, or user input s to define tissue type boundaries . Both hard callus and cortical bone are treated using the same piecewise material law and do not need to be separately segmented. This approach has key advantages for preserving objectivity because detecting the boundaries between old cortical bone and the most -mature hard callus is very challe nging using density -based thresholding (12). In ano ther study, we developed an image processing technique for assessing remodeling at the cortical -callus boundary in sheep. We found that during early healing, the cortical wall density decreases to create a smooth density gradient between the old and new bone (26). This evidence supports our decision to treat cortical bone and hard callus using a single material assignment la w and also highlights the potential challenges of pursuing a n alternative tri-zone material formulation for cortical bone, hard callus, and soft callus. Development of a region -specific scaling law for woven bone that is distinct from cortical bone would r equire density -based segmentation that can differentiate between old and new bone. In our experience with sheep, remodeling obfuscates the cortical -callus boundary such that cortical segmentation cannot be performed exclusively through density -based method s and may require user inputs that are not repeatable. Considering all the results together, the data provide support for the hypothesis that virtual mechanical testing with a dual -zone material assignment law that captures both the hard 12 and soft nature o f callus would out -perform single -zone material modeling when validated to physical biomechanical data. Good absolute agreement was achieved between all virtual and physical tests in Training and Testing groups (RMSEs lower than datum SD) . Correlations between the virtual and physical test results were slightly lower in the Testing group , but this was expected because the optimization was performed on the Training group . The dual -zone material model was still superior to the single -zone model at predicting biomechanical tests and it met the defined criteria for success for this material model validation. Overall, this suggests that when implemented with a dual -zone material model that captures the soft/hard heterogeneity of tissues within the callus , virtual mechanical testing is a reliable surrogate for physical biomechanical testing of healing ovine long bones. While the results of this study clearly show that modeling the dual soft -hard composition of the callus is necessary to achieve good absolute agree ment with physical validation data, the data do not provide a strong statistical justification to guide the choice of soft callus modulus , 𝐸𝑠𝑐. A literature review led to uncertainty regarding the true mechanical properties of fibrous and cartilaginous constituents of callus (27–33), so a three orders -of-magnitude variation in soft callus modulus value s was t ested (𝐸𝑆𝑐 = 0.5, 5, 50, or 500 MPa) . After optimization of the density cutoff for applying the soft modulus, we found that the dual-zone models with 𝐸𝑠𝑐 ≥ 5 MPa met our criteria for success based on quantitative analysis of torsional rigidity. This finding is both a strength and a weakness of our approach. It suggests that the bony architecture within the callus contributes most of the structural integrity and that modeling is relatively insensitive to the choice of soft tissue modulus when the objective is to assess whole -bone rigidity. However, defining success this way comes at a price, which is uncertainty regarding the localization of strains inside the callus. Considering both the quantitative analysis of 𝑅2 and RMSE and the 13 qualitative comparisons of strain maps and histology, we recommend implementation of the dual-zone material model with a soft callus modulus, 𝐸𝑠𝑐 = 50 MPa, for radiodensities below 665 mgHA/cm3. Further optimization of 𝐸𝑠𝑐 to achieve quantitative agreement with physical tests based on strain would require data that was not available in our mechanical tests , such as the use of image -based strain mapping during mechanical testing of sectioned samples (29). Unfortunately , these tests are destructive and have yet to be tried with a material -model optimization approach such as we have described, but a convergence of these methods merits further investigation. This work is not without limitations. CT scanning without contrast enhancement does not provide the ability to differentiate between soft tissues based on their radiodensity, so we were not able to use imaging to separ ately define regions of cartilaginous and fibrous soft tissues. We have also employed a linear elastic material model for all tissues in the dual -zone model. The use of linear elastic models for bone is widespread common practice in finite -element modeling of bone , but the soft tissues present in early fracture healing can have viscoelastic or poroelastic behavior (34–37). The mechanical test ing data for this study used a constant loading rate (5°/min) and produced linear torque -angle response curves (12). This observation lends support to the theory that the elastic bone contribution within the callus dominates the overall structural mechanics, at least at this stage of healing. Further research with non -destructive testing at multiple loading rates would be required to explore the suitability of piecewise nonlinear material models for representing these tissues. Finally, the translational potential of these findings must be considered in the context of our choice of animal model. O ur use of ovine data to develop the dual -zone material model is consistent with the widespread preference for sheep in osteotomy fracture models (38) due to 14 their similar bone healin g rate and remodeling potential in comparison with humans (39–41). However, sheep are not a perfect mechanical model for human bone. Ovine cortical bone is plexiform with few Haversian canals and its microstructure is more akin to woven bone than to human cortical bone (39). This may have contributed to our success in using a single material scaling law to assign mechanical properties across the range of radiodensities including hard callus and old corti cal bone. Hard or mature callus in humans may be more mechanically distinct from cortical bone than it is in sheep, but obtaining in vivo load-response data in healing human bones remains challenging, and this limitation may be difficult to overcome. In conclusion, this work shows how careful modeling of material inhomogeneity within the callus region is required to accurately replicate the measured torsional rigidity of osteotomized ovine tibiae using virtual mechanical tests . Use of a single material as signment law developed for intact cortical bone overpredicts the rigidity of healing bones because it treats all tissue within the healing zone as bone of varying density. In contrast, a dual -zone approach that differentiates between soft and hard callus a chieves strong correlations between virtual and physical test results with good absolute agreement . This material modeling method also reveals that much of the structural integrity of the healing bone is conferred from an internal architecture of mineraliz ed hard callus with interstitial soft tissue. Overall, the results demonstrate the translational potential of image -based virtual mechanical testing as a reliable, non -invasive assessment of bone healing. Materials & Methods Animal Study Information 15 Forty-four adult female Swiss alpine sheep (2 -3 years old, weighing 59 – 87 kg) were part of multiple previously completed research studies with three different tibial osteotomy models stabilized by medial plating (Fig. 5A). Detailed animal information for these studies has been previously reported (12). In brief, there were a total of 33 animals that had both postmortem biomechanical testing data and CT scans available for analysis. These animals comprised three different osteosynthesis mo dels, two of which were normal healing scenarios with 3 -mm osteotomies and one was a delayed healing scenario with a 17 -mm gap and autograft. 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 ) (12). Imaging, Mechanical Testing, & Scan Processing Micro -Computed Tomography (µCT) scans were performed after animal sacrifice using an XtremeCT II Micro‐CT scanner (Scanco Medical AG) with an X‐ray voltage of 68 kVp and X‐ ray current of 1470 µA. The resulting scans had an isotropic resolution of 60.7 µm. To convert from Hounsfield Units (HU) to calibrated radiodensity (𝜌𝑄𝐶𝑇, mgHA/cm3), a phantom was scanned with the same settings (Scanco KP70 phantom, QRM). All samples were tested in a custom -made fixture using an Instron E10000 electrodynamic testing machine. Axial loading and torque were measured with a calibrated load cell (±10 kN/ ± 100 Nm). The periosteum was stripped at the end sof each tibia, which wer e subsequently mounted into the custom fixture using Beracryl embedding medium. The periosteum was not stripped in the callus region . The diaphysis of each tibia was kept moist using saline‐soaked gauze. An axial preload of 5 N was maintained throughout the test. An internal rotation of 5° per minute was then applied. Biomechanical torsional rigidity was calculated as a linear regression 16 of the loading curve between 6 and 10 N-m multiplied by the distance between the proximal and distal surfaces of the Beracryl pots (the gauge length). The CT scans were processed using Materialize Mimics (21.0, Plymouth, MI ) following a workflow developed for a previous study (10). All µCT scans were first down -sampled to an isotropic resolution of 400 µm , which is comparable to clinical -resol ution scanning . Cortical bone and callus regions were segmented to create anatomically representative virtual models of the biological samples. Isometric 3D surface views of all 33 osteotomized models are provided in the SI Appendix , Fig. S4 . After segment ation , the masks were united into a single surface model in preparation for volumetric discretization. To produce cohesive surfaces, the united models were wrapped with a gap closing distance of 1mm and a smallest detail of 0.5 mm. Finite Element Modeling & Virtual Torsional Rigidity Testing After scan processing, t he finite element models were exported to 3-Matic (Mimics Innovation Suite ) and a quadratic tetrahedral (Tet -10) mesh was applied to each model with a maximum surface edge length of 1 mm and a maximum interior edge length of 1 mm. The meshing procedure and a supporting mesh convergence study were previously described (10). Elementwi se material properties were applied to the finite element meshes based on underlying scan radiodensity data using a piecewise material assignment law further outlined within the Material Model Definition section below. Virtual torsional rigidity (VTR) test ing was performed in ANSYS (2020 R2, Canonsburg, PA ). The distal end of each tibia was fixed in translation and rotation and the proximal end was rotated by one degree about the long axis of the bone, leaving all other degrees of freedom unfixed . For each model, virtual torsional rigidity was calculated as follows: 17 𝑉𝑇𝑅 =𝑀𝐿 𝜙 [2] where 𝑀 is the calculated moment reaction , 𝐿 is the working length of the test segment , and 𝜙 is the applied angle of twist . Callus Material Model Definition Fracture c allus is a heterogeneous, functionally graded composite structure comprised of a number of different tissue types , each of which have unique material properties (42). The premise of our dual-zone material model is that the tissues around the fracture line can be mechanically represented as having two distinct constitutive components: hard callus (bony regions of varying density ranging from newly mineralized woven bone up to cortical bone) and soft callus ( low-density fibrous and cartilaginous interstitial tissues). To accurately capture the structural mechanics of callus , it is necessary to treat these regions differently in the material model. To accomplished this, we defined a piecewise (dual -zone) material formulation that used a to-be-determined radiodensity threshold cutoff to differentiate between hard callus and soft callus. Above the cutoff, bone of var ying density was represented using a speci es-specific density -modulus relation. Previously, we used an optimization technique to derive a scaling equation for the density -dependent Young’s modulus of ovine cortical bone and achieve d close correspondence be tween physical and virtual torsion tests of intact ovine tibiae (12): 𝐸= 10225 ×(𝜌𝑄𝐶𝑇)1.0 [3] where E is Young’s modulus [MPa] and 𝜌𝑄𝐶𝑇 is the phantom -calibrated radiodensity [mgHA/cm3]. In the dual -zone material model, t his material assignment law for cortical bone was simply extended t o model the mechanical properties of all material (cortical bone and hard 18 callus) above the density cutoff threshold. A Poisso n ratio 0.3 was assumed for all elements (43, 44). To extend the single -zone material model of Equation 3 to a dual -zone piecewise application, two more parameters were required: a soft callus density cutoff (𝜌𝑐𝑢𝑡) and a soft callus modulus value (𝐸𝑠𝑐). The density cutoff value was used to specify the threshold below which elements were insufficiently mineralized to be treated as bone and should instead be modeled as soft tissue with modulus 𝐸𝑠𝑐. Below the cutoff, the soft tissues within the callus were treated as homogeneous with a modulus much lower than that of bone. After a literature search, a definitive value for 𝐸𝑠𝑐 could not be identified. However, upper and lower limits were identified based on previous work (27–33). The candidate soft callus modulus values were therefore chosen on a logarithmic range to be 𝐸𝑠𝑐= 0.5, 5, 50, and 500 MPa . These values were tested in combination with a range of 𝜌𝑐𝑢𝑡 density cutoffs v ia the procedure described in Material Model Optimization below to determine the piecewise material formulation best able to replicate in vitro biomechanical testing. Our density -delineated zonal approach to modeling callus was inspired by Mora -Macías et al. (2019), who used threshold -based segmentation to apply homogeneous mechanical properties within zones labelled as cortical bone, hard callus, and soft tissue (32). Their approach was adapted from earlier work by Shefelbine et al. (2005), who pos tulated a dual -zone piecewise approach for modeling material properties in murine fracture callus (33). We sought to combine the advantageous features of both methods (density -dependent mechanical properties of bone, density -based delineatio n between hard and soft callus) and seek a parameter set that would achieve close correspondence between physical and virtual mechanical tests of osteotomized ovine tibiae. 19 Material Model Optimization To perform parameter optimization for the dual -zone material model, t he N = 33 operated tibiae were randomly assigned to a Training group (N = 17) and a Testing group (N = 16) using permutations generated in MATLAB ( Fig. 5B). The optimization objective was to minimize the difference between calculated VTR and measured biomechanical torsional rigidity by varying 𝜌𝑐𝑢𝑡 and 𝐸𝑠𝑐 for all 17 Training group specimens simultaneously. The resulting response surface was evaluated based on root mean squared error (RMSE) between the virtual and physical torsional rigidities (Fig. 5C/D ). To assign material propert ies to an individual model for a candidate parameter set (combination of 𝜌𝑐𝑢𝑡 and 𝐸𝑆𝑐), first the raw material data for each element in the mesh was pre - processed using a custom function developed in MATLAB. In each finite element within the mesh, t he gray value density in Hounsfield Units was converted to calibrated radiodensity ( 𝜌𝑄𝐶𝑇, mgHA/cm3) using information from a calibration phantom scan. For each element, i f 𝜌𝑄𝐶𝑇 fell below the soft callus cutoff 𝜌𝑐𝑢𝑡, it was assig ned a soft callus modulus value (𝐸𝑠𝑐 = 0.5, 5, 50, or 500 MPa, as described above) . If 𝜌𝑄𝐶𝑇 for an element was above 𝜌𝑐𝑢𝑡, it was treated as bone and assigned a modulus using Equation 3 . The piecewise function can be represented as foll ows: 𝐸={𝐸𝑠𝑐 𝜌𝑄𝐶𝑇 <𝜌𝑐𝑢𝑡 10225 ×(𝜌𝑄𝐶𝑇)1.0𝜌𝑄𝐶𝑇 ≥𝜌𝑐𝑢𝑡 [4] To construct the response surface for optimization , testing was performed for each soft callus modulus value (0.5, 5, 50, and 500 MPa) with candidate val ues for the soft callus cutoff swept from 𝜌𝑐𝑢𝑡 = 0 to 1500 mgHA/cm3 in increments of 100 mgHA/cm3, resulting in 16 tested levels for the soft callus cutoff (Fig. 2B). At each candidate combination of 𝐸𝑠𝑐 and 𝜌𝑐𝑢𝑡, the 20 root-mean -square e rror (RMSE) between VTR and biomechanical torsional rigidity GJ was calculated based on all 17 Training group samples. Across the range of 𝜌𝑐𝑢𝑡 values , the approximate minimum RMSE was then identified for each soft callus modulus (0.5, 5, 50, and 500 MPa). The virtual torsion tests were then run on the Training group models at additional 𝜌𝑐𝑢𝑡 refinement points on either side of the approximate minima in increments of 5 mgHA/cm3 (Fig. 2C). After refinement, the RMSE minima were identified for e ach of the four candidate 𝐸𝑠𝑐 value s and the four corresponding optimized 𝜌𝑐𝑢𝑡 values were recorded ( Fig. 2D). Altogether, construction of the refined response surface for all models in the Training group required a total of 2,363 virtual mech anical tests. All simulations were run on nodes of the High -Performance Computing (HPC) server hosted on Lehigh University’s campus and required a total of approximately 5,000 CPU -hrs. To test the validity of the optimized piecewise material model, each o f the four optimized soft callus cutoffs were then applied to the N = 16 Testing group models. The RMSE for VTR versus GJ was then calculated at each of the optimized combinations of 𝐸𝑠𝑐 and 𝜌𝑐𝑢𝑡 and the structural simulations results were saved for performance comparison between each of these four states . Data Processing, Visualization, and Statistic s Methods Simulation post -processing was performed in MATLAB. Mesh data, material model data, and the mechanical testing results from ANSYS were queried using custom MATLAB functions. Contour slice plots and 3D visualizations of models and results were generated using Paraview (5.9.0 , Sandia National Laboratories, Kitware Inc, Los Alamos National Laboratory) . 21 Descriptive statistics were generated using Microsoft Excel. Additional statistical analyses were generated using IBM SPSS Statistics 27 (Armonk, NY ). Pearson’s correlations were used to determine the strength of associatio n between virtual and biomechanical torsional rigidity datasets. A correlation coefficient 𝑅 ≥ 0.8 ( 𝑅2 ≥ 0.64) was defined as strong and 𝑅 ≥ 0.6 (𝑅2 ≥ 0.36) was defined as moderate (45). The criteria for evaluating “good” absolute agreement between correlated measures was the requirement that r oot mean squared error (RMSE) be less than the standard deviation of the datum set. One-way repeated -measures ANOVA was performed on both the Training and Testing groups to determine if there were statistically significant differences between the measured GJ from biomechanical testing and predicted VTR in each of the four optimized material model s (combination s of 𝐸𝑠𝑐 and 𝜌𝑐𝑢𝑡). Post hoc analysis with a Bonferroni adjustment was performed on pairwise comparisons to identify statistically significant differences . All values reported are averages and standard deviations unless otherwise stated. The statist ical significance limit was p = 0.05. Acknowledgments This material is based upon work supported by the National Science Foundation (NSF) under Grant No. CMMI‐1943287. The preclinical studies from which this data was obtained were 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 studies data for our analyses. Porti ons of this research were conducted on Lehigh University's Research Computing infrastructure partially supported by NSF Award 2019035. PS was supported in part by an NSF / 22 Institute for International Education (IIE) Graduate International Research Experien ce (GIRE) Fellowship. The authors declare no competing interests directly relevant to this work . BvR was the head of MSRU at the time of the animal experiment s and had oversight of these and other studies funded by Johnson & Johnson Family of Companies. HLD has stock or stock options, receives research or institutional support as a principal investigator, and receives financial or material support from OrthoXel. Additionally, HLD serves on a committee in the Orthopedic Research Society . 23