(150 -250 words : currently 2 22) 17 Bone healing has been traditionally descri bed as a four -phase process : inflammatory response, soft callus 18 formation, hard callus development, and remodeling. The remodeling phase has been largely neglected in most 19 numerical mechanoregulation models of fracture repair in favor of capturing early healing using a pre -defined 20 callus domain. However, in vivo evidence suggests that remodeling occurs concurrently with repair and causes 21 changes in cortical bone adjacent to callus tha t are typically neglected in numerical models of bone healing. The 22 objective of this study was to use image processing techniques to quantify this early -stage remodeling in ovine 23 osteotom ies. To accomplish this, we developed a numerical method for radiodensity profilometry with 24 optimization -based curve fitting to mathematically model the bone density gradients in the radial direction across 25 the cortical wall and callus. After assessing data from 26 sheep, we defined a dimensionless density fitting function 26 that revealed significant remodeling occur ring in the cortical wall adjacent to callus during early healing, a 23 % 27 average reduction in density compared to intact. This fitting function is robust for modeling radial density 28 gradients in both intact bone and fracture repair scenarios and can capture a wide variety of the healing responses . 29 The fitting function can also be scaled easily for comparison to numerical model predictions and may be useful for 30 validating future mechanoregulatory models of coupled fracture repair and remodeling. 31 32 Keywords: Rem odeling, Bone fracture healing , Micro -computed tomography , Tibia , Ovine osteotomy 33 Acknowledgments 34 The preclinical studies from which this data was obtained were funded by the Johnson & Johnson Family 35 of Companies. The authors wish to thank Beat Lechmann (Johnson & Johnson Family of Companies) for agreeing 36 to grant access to the ovine study data for our analyses. This material is based upon work supported by the 37 National Science Foundation under Grant No. 1943287. Any opinions, findings, and conclusions or 38 recommendations expressed in this material are those of the author(s) and do not necessarily reflec t the views of 39 the Johnson & Johnson Family of Companies or the National Science Foundation. 40 Conflicts of interest 41 BvR is a paid consultant and receives research support as a principal investigator from Johnson & 42 Johnson Family of Companies. The other auth ors report no conflicts relevant to this work . 43 Introduction 44 The fracture healing proc ess is commonly described as a four -phase concept: inflammation phase, soft 45 callus formation phase , hard callus formation phase , and bone remodeling phase (Claes et al. 2012; Pivonka and 46 Dunstan 2012; Einhorn and Gerstenfeld 2015; Ghiasi et al. 2017) . Numerical simulations have been used to predict 47 mechanobiological processes in the regions of fracture healing. Numerous finite -element (FE) models have been 48 developed to analyze local mechanical responses to varying loading conditions (Ament and Hofer 2000; Claes et 49 al. 2002; Lacroix and Prendergast 2002; Shefelbine et al. 2005; Simon et al. 2011; Steiner et al. 2013, 2014; Wilson 50 et al. 2017; Wang and Yang 2018; Pietsch et al. 2018; Ghiasi et al. 2019; Engelhardt et al. 2021) . Such models have 51 been used to study different factors that impact the healing process , for example the mechanics of the fixation 52 hardware, or outcomes following inverse dynamization (Wehner et al. 2014; Wilson et al. 2017) . Some studies 53 have also emphasize d that the initial phase of healing is critical for modeling of the healing process (Ghiasi et al. 54 2019) . Recently, we applied a mechanoregulatory model of fracture healing to show that the rate of healing 55 depends on the geometry of the fracture and that all fracture types experience a small healing delay with torsional 56 instability (Ren and Dailey 2020) . We also show ed that modifying a strain -based mechanoregulatory model of 57 fracture healing to include a geometric proximity function can successfully restrict callus formation to the region 58 near the fracture line and achieve simulation domain independenc e (Schwarzenberg et al. 2021) . 59 These techniques hold promise to predict outcomes under differ ent mechanical or biological conditions 60 and in response to new treatment strategies , but they all share a key limitation : the cortical wall is assigned static 61 mechanical properties and included in the mechanical loading simulations, but it is excluded from the 62 mechanoregulatory healing algorithms. T hese models neglect remodeling activity during fracture healing because 63 this last stage of fracture healing is often described conceptually as having limited overlap with the repair phase . 64 However, in vivo data suggests that remodeling of the cortical wall occurs concurrently with soft callus 65 development and mineralization. For example, p revious reports generated from nanoindentation testing have 66 provided highly localized evidence of decreasing bone density in the radial direction through the cortical bone 67 toward the callus (Preininger et al. 2011) . This evidence indicates the remodeling phase of secondary bone 68 fract ure healing is associated with structural changes in both the callus and the adjacent cortical bone throughout 69 the healing process and that these effects should be included in early -stage fracture repair models . 70 The key barrier to developing new mechanoregulatory models that can consider coupled repair and 71 remodeling is the lack of available validation data. The validation strategy used for most recent 72 mechanoregulation models has focused on the global structural mechanics of the limb and qualitative comparison 73 of model predictions to histological analysis of tissue in the callus . For example , the assessment of a model ’s 74 performance in predicting the spatial and temporal development of tissues in the callus is usually based on only 75 the general trends of mechanical data (e.g., torsi onal rigid ity) (Ren and Dailey 2020; Schwarzenberg et al. 2021) 76 or temporal data (e.g., bridgin g time, interfragmentary motion ) (Lacroix and Prendergast 2002; Isaksson et al. 77 2006; Wehner et al. 2014; Pietsch et al. 2018) . Ideally, t he localized tissue composition and distribution should 78 also be investigated to assess whether tissue growth and differentiation criteria within the algorithms are 79 functioning correctly . The historical gold standard for assessing mechanical properties within tissue samples has 80 been nanoindentation testing. This method has been used to probe the micro -scale properties of bone and callus 81 and identify gradients in porosity and density (Manjubala et al. 2009; Preininger et al. 2011; Mora -mací as et al. 82 2017) . However, these physical measurement techniques involve time -consuming destructive postmortem testing 83 with specialized equipment and the resulting data cannot be obtained in vivo , or for large regions of tissue, or in 84 large numbers of samples. Regrettably, these limitations mean that generalizable quantitative data describing 85 remodeling across the cortical wall into callus does not currently exist. 86 Accordingly, the purpose of this study was to develop a numerical technique for noninvasively quantif ying 87 structural remodeling adaptations at the bone -callus boundary in sheep by investigating the radiodensity data 88 from micro -computed tomography (µCT) scans . We hypothesize d that measurable cortical remodeling happens 89 early in the healing period, in conjunction with bridging of the external callus in secondary fracture he aling. 90 91 92 Methods 93 Animal Specimen Information 94 Twenty -six adult female Swiss alpine sheep (2 -3 years old, weighing 59 – 87 kg) were part of two previously 95 completed research studies with two different tibial osteotomy models stabilized by medial plating. Grou p 1 96 consisted of data from twelve animals with a 17 mm tibia defect augmented with autografts and stabilized with a 97 13-hole stainless steel plate (broad straight veterinary 3.5 mm LCP, 172 mm in length, with 3.5 mm bicortical 98 screws; DePuy Synthes®). Sheep in group 1 were sacrificed 12 weeks after surgery . Group 2 consisted of data from 99 fourteen animals with a 3 mm gap tibia defect stabilized with a 12 -hole stainless steel plate (broad straight 100 veterinary 3.5 mm locking compression plate (LCP), 159 mm in len gth, with 3.5 mm bicortical screws; DePuy 101 Synthes®). Sheep in group 2 were sacrificed 9 weeks after surgery. All experiments were conducted at the 102 Musculoskeletal Research Unit in Zürich, Switzerland , according to the Swiss laws of animal protection and we lfare 103 and authorized by the local governmental veterinary authorities (License No. ZH 183/17) . 104 Micro -Computed Tomography (μCT) Scanning 105 After animal sacrifice, both intact and osteotomized tibiae were excised, stripped of soft tissue, and all 106 hardware was removed , taking care not to disrupt the callus. Samples were then wrapped in saline -soaked gauze, 107 and µCT scanned using an XtremeCT II Micro -CT scanner (Scanco Medical AG, Bruettisellen, Switzerland) with an 108 X-ray voltage of 68 kVp and X -ray c urrent of 1470 µA. The resulting scans had an isotropic resolution of 60.7 µm. 109 Image Segmentation 110 A region of interest (ROI) was selected at three different location s of each operated limb diaphyseal 111 segment scan : 10 mm from the most proximal aspect of the tibia , 10 mm from the most distal aspect of the tibia , 112 and just proximal to the fracture line at the level of maximum callus radius (Fig. 1). The proximal and distal ROI 113 were above and below the most proximal and distal screw holes . Sampling slices from the intact limb were also 114 taken at the corresponding level s. At each sampling location (region of interest), a total of 20 slices were selected 115 for analysis . 116 An automated segmentation algorithm was developed in MATLAB 2021a (The MathWorks, Inc., Natick, 117 Massachusetts, USA) to define the outer boundary of the callus as well as the area enclosed by the pericortical 118 boundary of the preexisting cortical bone . The detailed workflow can be seen in Fig. 2. Post-segmentation data 119 processing was carried out on all 2D tomogram s located within each ROI. 120 The process for s egmentation of cortical bone and callus was based in part on previously reported 121 methods (Buie et al. 2007; Morgan E, Mason Z, Chien K, Pfeiffer A, Barnes G, Einhorn T et al. 2009; Bissinger et al. 122 2017) . To complete this task, the algorithm required two initial estimat es of the segmentation thresholds . The 123 selection of these thresholds affect s the cortical bone and callus boundaries of the output mask and require s 124 subsequent refinement . First, the soft tissue region was exclude d using connected -component labeling . Noise was 125 then reduced by applying a median filter. Next , the original dataset was re-threshold ed and masked with the 126 detected boundaries of the bone -callus region. Dilation and erosion operations were then used to determine the 127 bone and callus boundary . The initial estimat es of the segmentation threshold s were then refined by minimiz ation 128 of spline energy to achieve smoothness of the detected tissue boundar ies for each sheep. A local ly convex 129 boundary criterion was also used to correct any boundary defects caused by small region s of lower-density cortical 130 bone. Finally, a mask of the callus , cortical bone , and non -bone regions was created by combining the output of 131 the previous steps. 132 The slice extraction and segmentation processes were repeated at all three level s of height in both the 133 intact and operated tibia e. No callus was identified in the intact limbs. 134 Radiodensity Data Sampling and Curve Fitting 135 The process for radiodensity sampling is illustrated in Fig. 3. In each 2D tomogram within the region of 136 interest, t he centroid of the cortical bone was calculated from the endocortical boundary (red boundary) in the 137 mask geometry. A polar array of sampling lines was created in 1-degree increments throughout the entire 360 - 138 degree angle range, with each line having its origin at the centroid and extending through the cortical wall and 139 callus. The radiodensity values along the profile lines were sampled in 60 µm increments in the radial direction , a 140 distance equal to the in -plane resolution of the scan . Pixel data was retrieved from the grayscale image using a 141 bicubic interpolation technique . Data quality was checked for each profile line and individual p rofile line s were 142 excluded from the data set if the callus was not present or too thin ( thickness smaller than 10mm , measured in 143 the radial direction along a profile sampling line ). Radial sampling was repeated using the same procedure for all 144 20 slices in each ROI and pooled to create a point cloud that defined the radiodensity profile line for each ROI in 145 each animal ( Fig. 1c ). 146 In this animal model, medial plating produce d asymmetric callus formation, as can be clearly seen in Fig. 147 3c. Cortical wall thickness was also not constant across all radial profile lines. To enable a combined curve -fitting 148 analysis of the profilometry samples from within each region of interest, a registration and scaling procedure was 149 defined to create a consistent dimensionless coordina te system for the cortical wall and for the callus, when 150 present. Each profile line was registered to set its linear coordinate system to zero at the bone -callus boundary , 151 or the pericortical boundary if intact. Negative linear coordinates corresponding to the cortical wall were 152 normalized to the local cortex thickness of each profile line (Fig. 3c/d), such that 𝑟 = -1 was the endocortical 153 boundary and 𝑟 = 0 was the pericortical boundary . For the positive linear coordinates, the normalization 154 procedure de pend ed on whether callus was present. For the intact limb s, positive coordinates were normalized 155 to the local cortex thickness (Fig. 5b ). For the operated limb s with callus, positive coordinates were normalized to 156 the doubled local callus thickness , producing 𝑟 = 0.5 at the outer boundary of the callus mask and ensuring data 157 sampling beyond the callus into the soft tissue zone ( Fig. 5c). The resulting average number of sampling points 158 per region of interest was 170,340. 159 To enable compariso ns between the treated and untreated limbs of all animals, the sampled r adiodensity 160 values were also normalized to the median cortical bone density of the intact diaphysis of each animal. Within 161 each animal, all local radiodensity values sampled from the profile lines were normalized to this baseline median 162 cortical bone density , regardless of limb (intact or operated) and ROI location (proximal, midshaft, and distal) . 163 Curve fitting techniques w ere then used to mathematically model the collective radial density profile using 164 the sampling points collected at each ROI. To develop an appropriate fitting function for measuring local structural 165 changes, we started with the simplest case , an intact cortical wall fit with a double sigmoid funct ion: 166 𝑓(𝑟)=1 1+𝑒−𝑎1(𝑟−𝑏1)∗1 1+𝑒−𝑎2(𝑟−𝑏2) (1) where in general for a sigmoid function, the parameters 𝑎 sets the width of the transition zone from a function 167 value of 0 to a function value of 1, and the parameter 𝑏 sets the center of this transition zone. 168 Next , a logical linear term and a constant were added to the double sigmoid function to account for the 169 possibility of a nonzero density gradient in the cortical wall ( 𝑟<0) and nonzero background radiodensity : 170 𝑓(𝑟)=1 1+𝑒−𝑎1(𝑟−𝑏1)∗1 1+𝑒−𝑎2(𝑟−𝑏2)+𝑑+{𝑐𝑟 𝑟<0 0 𝑟≥0 (2) where 𝑐 is the linear slope applied only within the cortical wall and 𝑑 is the background density constant. 171 To adapt this method for profilometry samples with callus, an additional sigmoid term was added : 172 𝑓(𝑟)={ 1 1+𝑒−𝑎1(𝑟−𝑏1)∗(𝑔1+(1−𝑔1) 1+𝑒−𝑎2(𝑟−𝑏2))+𝑐𝑟+𝑑1 ,𝑟<0 𝑔2 1+𝑒−𝑎3(𝑟−𝑏3)+𝑑2,𝑟≥0 (3) where the constants 𝑔1 and 𝑔2 in the sigmoid terms allow for non -unity in the upper sigmoid plateau and the 173 constants 𝑑1 and 𝑑2 allow for nonzero background radiodensity on both the internal and external bounds of the 174 fitting function. In this piecewise definition, the cortical fun ction (𝑟<0) is a generalized double logistic function, 175 and the callus function is a single logistic function. 176 Parameter values for constants in Eq. 2 -3 were fit using the nonlinear least -squares method via the 177 MATLAB Curve Fitting Toolbox (R2021a). I n addi tion to this standard curve -fitting technique, a multi -objective 178 optimization method was introduced to fit the parameters of the function defined in Eq. 3 with additional 179 constraints requiring position and slope continuity of the piecewise fit function at 𝑟=0 as described in Eq. 4 . This 180 multi -objective optimization method can be formulated to find the coefficient s of 𝑓(𝑟), minimizing the residual s 181 (sum of the square s), of all the data point s: 182 183 184 185 186 𝑚𝑖𝑛 𝑎,𝑏,𝑔,𝑑∑(𝑦𝑖−𝑓𝑟<0(𝑟𝑖))2𝑚 𝑖=1 𝑚𝑖𝑛𝑎,𝑏,𝑔,𝑑,𝑐∑(𝑦𝑘−𝑓𝑟>0(𝑟𝑘))2𝑛 𝑘=1 𝑠.𝑡.𝑓𝑟<0(0) −𝑓𝑟>0(0)=0 𝑓𝑟<0′(0)− 𝑓𝑟>0′(0)=0 (4) 187 where 𝑟𝑖 and 𝑟𝑘 are the set of data points in the cortical wall or external to the cortical wall, 𝑛 and 𝑚 are the total 188 number of data points of the cortical wall or the callus. The minimization problem was solved using the MATLAB 189 Optimization Toolbox ( R2021a ). After the fitting parameters were obtained for each sheep, the resulting fit 190 function, 𝑓(𝑟), was evaluated to calculate the slope at the bone -callus interface to quantify the extent to which 191 remodeling had produced a smooth functionally graded interface at the bone -callus boundary. For intact bone 192 not bounded by callus, an additional custom MATL AB script was also used to calcul ate the average density value 193 in the cortical wall plateau , which was reported as the m ean cortical density of each fitted curve. 194 Statistical Analysis 195 All statistical analyses were performed using SPSS Statistics 27 (IBM Corp; Armonk, New York, USA). A 196 value of P = 0.05 was used as the threshold for statistical significance. A two -tailed paired-sample t-test was used 197 to analy ze if the fitting function co nstants and the evaluated slope of the fitting function at the pericortical 198 boundary were significantly different between the o perated and contralateral limbs. 199 Results 200 Choice of Curve Fitting Function 201 Several alternative methods were investigated and critically evaluated with respect to their ability to 202 capture the essential features of the bone density profile curves. The fitting results of the double sigmoid function 203 (Eq. 2) are shown in Fig. 4a. This function provided a reasonable a pproximation for the intact cortical wall, but it 204 was not able to capture the increasing slope of the density in the radial direction for intact cortical bone and it 205 missed both the cortical and callus plateaus in the data for the operated tibia . Adding a logical linear term and a 206 constant produced a good fit with 𝑅2≥0.93 for the double sigmoid function (Eq. 3) with cortical bone , but this 207 fitting function fail ed to capture the remodeling behavior evidenced by a decrease in density in the cortical wall 208 at the bone callus boundary ( Fig. 4 b). The p iecewise logistic function (Eq. 4) was then used to get a better fit when 209 callus was present , this function provide d a better fitting quality ( 𝑅2≥0.95), but when fit using a nonlinear least 210 squares minimization, the function was not continuous at the bone callus interface (Fig. 4c ). 211 For both the intact and operated tibia e, the best fitting method was to use the multi -objective 212 optimization method (Fig. 4d ) to identify the parameters of the fitting function shown in Eq. 4. This produced 213 continuous fit functions with goodness of fit statistics show ing an overall 𝑅2≥0.87 for all intact tibiae and 214 operated tibiae . In all samples, all six fitting parameters had 95% confidence intervals that did not include zero. 215 The fitting parameters for all 26 animals at each of the six sampled ROIs can be seen in the supplementary digital 216 content . 217 The performance of the final fitting function (Eq. 4 with multi -objective optimization) is shown for a 218 different animal in Fig. 5. This single function successfully captured several observed characteristics of both the 219 intact and operated limb. Specifically, in the intact tibia, it captured the observ ed increase in radiodensity in the 220 radial direction from the endosteal to the pericortical boundary (Fig. 5a). For the operated limb, the fit curve 221 captured the overall decrease in radiosity in the cortical wall and the localized remodeling at the pericort ical 222 boundary bounded by callus ( Fig. 5b). 223 Each paramete r of the fitting function was included to capture the observable features of the profilometry 224 data in different region s of the fitted curve and at transitions between regions. The resulting constants can then 225 be associated with physical features related to the bone -callus radial density function . The magnitude s of 226 parameters 𝑎1 and 𝑎2 define the width of the region of density change at the perioste al and endocortical 227 boundary s of cortical bone . The parameter s 𝑏1 and 𝑏2 define the location of the center of the transition of the 228 perioste al and endocortical boundary s of cortical bone. Parameter 𝑎3 define s the width of the region of density 229 change at the outer boundary of the ca llus, while parameter 𝑏3 sets the location of the center of the transition 230 between callus and the surrounding soft tissue . The linear term c in the cortical bone zone ( 𝑟 ≤ 0) defines the 231 increasing density of cortical bone from the endosteal to pericortical boundary . Finally, the c onstant 𝑑1 capture s 232 the average density of cortical bone and the term 𝑑2 compensate s for the noise of the data outside the callus. 233 Finally, 𝑔1and 𝑔2 capture the cortical density at the pericortical boundary of the injured limb . 234 235 236 Cortical Tissue Alterations in Operated Limbs 237 In addition to profilometry sampling at the level of the fracture, we also analy zed cortical bone at the 238 proximal and distal ends of the diaphysis of the intact and operated limbs. The resulting profilometry curves of all 239 animals combined are shown in Fig. 6. Comparing the operated group to their matched intact tibiae at both the 240 dista l and proximal ends of the bones , there were several changes in the cortical walls of the operated tibiae, both 241 at the level of the fracture and at the proximal and distal ends of the bone . Table 1 shows the comparisons 242 between the operated and intact limb s based on two summary parameters derived from the fitting functions — 243 density slope at the pericortical boundary and mean cortical wall density. Each parameter was evaluated at the 244 proximal , midshaft, and distal ROIs in both the intact and osteotomized limb s. Both parameters were significantly 245 different between operated and intact for all ROIs, although the cortical density decreases in the operated limb 246 were 10x larger at the midshaft than at the proximal and distal ends of the bone. Considering all animals, sampling 247 at the proximal and distal ends of the bone showed that t he mean cortical density was 1.3% and 2.2% lower 248 respectively in the operated limbs compared to intact (p < 0.001). 249 All animals also showed evidence of remodeling in the cortical wall bounded by callus (Fig. 6d ). The 250 preexisting increase in radiodensity toward th e pericortical boundary that can be seen in the intact bones ( Fig. 6a ) 251 disappeared in the osteotomized bones ( Fig. 6 c). At the pericortical boundary bounded by callus, these changes 252 corresponded to a 23.4 % drop on average from the peak density of the native cortical wall. Overall, comparing 253 the intact and operated groups at the midshaft , there were statistically significant differences in all obtained 254 radiodensity profile fitting function parameters except for 𝑔1 (p ≤ 0001; see supplementary Table S1). The slope 255 of the fitted density profile curve at the pericortical boundary was also significantly reduced in operated limbs 256 compared to intact at all three ROIs (p ≤ 0.007; see Table S1). 257 Discussion 258 This study presents a scalable mathematical method for non -invasive quantification of tissue 259 mineralization density of the cortical tissue matrix and callus within a n ovine tibia l shaft. The results show that 260 remodeling alterations occur not only in the callus tissue due to mineralization but also in the cortical tissue 261 adjacent to the repair zone thro ughout the healing process . The findings are consistent with expectations from 262 destructive physical testing and confirm the presence of a decrea sing gradient in mineral density of the cortical 263 bone at the vicinity of a fracture during secondary healing. We previously showed that at 9 weeks, most of these 264 osteotomies had achieved a structural rigidity comparable to intact tibiae (Schwarzenberg, 202 0). This 265 profilometry curve fitting data shows that in sheep, significant adaptive changes have occurred in cortical bone 266 bounded by callus , concurrent with structural bridging as early as 9 or 12 weeks post -op. It should be noted that 267 the remodeling of co rtical bone was measurable not only near the fracture site, but also to a lesser extent at the 268 proximal and distal end of the operated limbs. The methods described herein could be use d to measur e 269 remodeling as a temporal in vivo outcome measure from imaging only. These measures could be useful as an 270 indicator of healing speed by establishing the status of a fracture on the continuum of mechanical and structural 271 adaptations that arise during healing. 272 The intended purpose of this investigation was to provide a useful set of validation data for the future 273 development of mechanoregulatory models that include the cortical wall and address coupled remodeling at the 274 cortical -callus boundary. The results confirm that significant changes do occur on the cortical side of the cortical - 275 callus boundary and that not addressing these changes is a critical limitation in the current comput ational models 276 of bone fracture healing . These models usually assume a time -invariant value for cortical bone stiffness, which 277 does not reflect the in vivo reality of a 23% drop in cortical radiodensity at the cortical -callus boundary . The 278 normalized generalizable fitting function we reported could be used to assess whether mechanoregulatory models 279 that include th e cortical zone are capturing the appropriate magnitude of changes, and whether the distribution 280 of predicted material properties at the boundary and within the callus region are valid . 281 Investigators who are interested in applying our normalized, dimensio nless remodeling curve to validate 282 their dimensional mechanoregulation models may wish to use the following scaling procedure. First, observe that 283 the fitting curve cannot be used directly because the function has been normalized in both the radius coordin ates 284 (𝑟) and radiodensity values ( 𝑦 axis). To dimensionalize the radial coordinates, recall that the radial coordinates 285 were normalized to the cortical wall thickness for 𝑟<0 and to the doubled callus thickness for 𝑟 > 0. To scale this 286 function appropriately in 𝑟, identify the cortical w all thickness, 𝑡, and callus thickness, 𝑡𝑐. The y coordinates in our 287 fit curve were normalized to the median radiodensity [HU] of the intact limb for each animal. For the N = 26 sheep 288 analysed in this study, the resulting average intact tibia radiodensi ty was 3164 ± 10 HU. The apparent radiodensity 289 [HU] must be scaled to bone mineral density by means of a linear conversion derived from a phantom scan. This 290 density value can then be converted to a Young’s modulus using a species -specific scaling law. For ovine cortical 291 bone, we recommend the equation derived by (Schwarzenberg et al. 2021) , which produces an average Young’s 292 modulus, 𝐸𝑖𝑛𝑡𝑎𝑐𝑡 = 12.9 GPa for intact bone. After selecting 𝐸𝑖𝑛𝑡𝑎𝑐𝑡 and matching 𝑡 and 𝑡𝑐 to the chosen 293 computational model, the following scaled remodeling function can be used to validate the model predictions: 294 𝐸(𝑟)=𝐸𝑖𝑛𝑡𝑎𝑐𝑡∗ { (1 1+𝑒−𝑎1(𝑟 𝑡−𝑏1)∗(𝑔1+(1−𝑔1) 1+𝑒−𝑎2(𝑟 𝑡−𝑏2))+𝑟𝑐 𝑡+𝑑1) ,𝑟<0 (𝑔2 1+𝑒−𝑎3(𝑟 𝑡𝑐−𝑏3)+𝑑2),𝑟≥0 (5) Where the fitting constants ( 𝑎1,𝑏1, 𝑔1, 𝑎2, 𝑏2, 𝑐, 𝑑1, 𝑎3, 𝑏3, 𝑔2, and 𝑑2) are drawn from Table S1. This 295 dimensionalized remodeling curve will have radial coordinates 𝑟 that are zeroed at the cortical -callus boundary, 296 but are dimensionally scaled to match the physical size of the geometry in the computational model. 297 This study is not without limitations . First, the sheep imag ing data were only available at the time of 298 sacrifice . As a result, no information about the temporal variations of local mechanical properties of the cortical 299 tissue is currently available for this study . However, due to the inclusion of both large -gap and small -gap healing 300 models, the selected sheep groups represent a wide range of healing responses and callus sizes, suggesting that 301 the numerical method s and resulting fitting function we have proposed are robust . Another key limitation is that 302 we applied the segmentation algorithm to evaluat e only a small fraction of the 2 D tomograms that comprise the 303 entire cortical wall and focused on gradients in the radial direction . Remodeling effects in the cortical wall also 304 have an axial dependency on proximity to the fracture line, which has not been explored here. Additional work 305 would be required to achieve profilometry analysis of tissue density profile s in a more general 3D sense. Finally, 306 it is important to note that the resolution of the medical images themselves could also influence the curve fitting 307 quality. In this study, we have used high -resolution of micro -CT scans that were obtained post -mortem. Applying 308 these same techniques to in vivo computed tomography wou ld likely require larger voxel sizes and may not 309 provide enough data points for curve fitting purpose in the area of transition especially at the cortical -callus 310 boundary. We recommend that future studies explore whether larger voxels are inferior for dete ction of 311 remodeling effects in the cortical wall. 312 Conclusions 313 This study show ed that structural remodeling at the bone -callus interface can be detected in sheep using 314 radiodensity profilometry from micro -CT scans. It also show ed that the density of cortical bone adjacent to callus , 315 which correlates with its elastic properties, is not constant during the early healing period and that this region 316 represents an integral part of the process of coupled fracture repair and remodeling. The scalable function 317 presented herein for mathematically representing density gradients in the cortical wall and callus may be useful 318 for validating future mechanoregulatory models of fracture repair . 319 Figure Captions 320 Figure 1: Fig. 1 – (a) Coronal slice view from a micro -CT scan of an ovine tibial osteotomy at 9 weeks post -op. In 321 this study, radiodensity gradients were sampled in the direction indicated by the red arrow. (b) Three regions of 322 interest (ROI) were identified for sampling in each intact and each operated limb: proximal, distal, and at the level 323 of the fracture. (c) Each ROI consisted of 20 axial slices that were analyzed to measure radial radiodensity profiles 324 in the cortical wall and callus region. 325 326 Figure 2: Regions of cortical bone and cal lus were identified using an automated segmentation method developed 327 in MATLAB 2021a. This method requires two global threshold inputs, which are used to initialize the splines that 328 define the endosteal, peri cortical , and callus boundaries, and subsequentl y refined as shown. 329 330 Figure 3: (a) A polar array of radial profile lines was created from the cortical centroid through the cortical wall 331 and callus, if present, and radiodensity values were sampled from the underlying scan data. (b) For intact bones, 332 each profile line was registered to set its linear coordinate system to zero at the pericortical boundary and 333 sampling points were extended up to one cortical thickness into the surrounding soft tissue. (c) For samples with 334 callus, profile lines were registere d to zero at the bone -callus boundary and sampling points were extended to 335 double the local callus thickness. (d) Resulting compiled profilometry sampling data for one animal at the level of 336 the callus in the osteotomized limb. All radiodensity values were normalized to the median cortical bone density 337 of the intact tibia for that animal. 338 Figure 4: Profilometry data sampled at the level of the fracture for intact tibiae (black, left) and operated tibiae 339 (red, right). Fit curves from the three tested fitting equations in panels (a) -(c), with all parameters obtained 340 following a least -squares minimization approach using the Curve Fitting Toolbox in MATLAB. Fit curves in panels 341 (c) and (d) both use Eq. 4, but additional continuity constraint at 𝑟 = 0 was achiev ed in (d) using multi -objective 342 optimization. 343 Figure 5: Profilometry sampling points and resulting fitted curves using the method of Fig. 4(d) in a different animal. 344 The radial coordinate registration procedure allowed consistent referencing of key physica l locations across 345 different animals: endocortical boundary , pericortical boundary , and callus outer boundary (if present). The 346 density normalization procedure eliminated baseline differences in bone properties between animals by scaling 347 all sampling data points to the median cortical bone density for the intact tibia in each animal. 348 349 Figure 6: Profilometry curves were fitted from data sampled at the same levels at the proximal end (a) or distal 350 end (b) end of intact tibiae (black curves) and operated tibiae (red curves). Profilometry curves were also fitted 351 from data sampled at the midshaft of all intact tibiae (c, black curves) and operated tibiae (d, red curves). 352 353 Table 1: Comparison of profilometry curve fitting parameters between operated and contralateral groups . 354 355 356 357 358