Subject-specific mechanical loading of bones: Improved healing of mouse vertebral defects using adaptive loading regimes

Methods to repair bone defects arising from trauma, resection, or disease, continue to be sought after. Cyclic mechanical loading is well established to influence bone (re)modelling activity, in which bone formation and resorption are correlated to micro-scale strain. Based on this, the application of mechanical stimulation across a bone defect could improve healing. However, if ignoring the mechanical integrity of defected bone, loading regimes have a high potential to either cause damage or be ineffective. This study explores real-time finite element (rtFE) methods that use three-dimensional structural analyses from micro-computed tomography images to estimate effective peak cyclic loads in a subject-specific and time-dependent manner. It demonstrates the concept in a cyclically loaded mouse caudal vertebral bone defect model. Using rtFE analysis combined with adaptive mechanical loading, mouse bone healing was significantly improved over non-loaded controls, with no incidence of vertebral fractures. Such rtFE-driven adaptive loading regimes demonstrated here could be relevant to clinical bone defect healing scenarios, where mechanical loading can become patient-specific and more efficacious. This is achieved by accounting for initial bone defect conditions and spatio-temporal healing, both being factors that are always unique to the patient.


Introduction 31
The management of critical-size bone defects continues to present surgical challenges. Trauma and 32 bone resection can lead to lengthy recovery times or amputation. The use of autografts is the current 33 gold standard, however, is quantity-limited and accounts for 20% of the complications 1 . Despite 34 advances in biomaterial development and understanding of signaling mechanisms, the search for 35 improved treatment methods of such bone defects continues. 36 37 One treatment method of interest is mechanical loading of bone. Mechanical interactions have a long-38 established relationship to bone physiology, leading to earlier concepts of micromotion during bone 39 defect healing 2 , and the significance of fracture instability in healing outcomes 3,4 . The effects of 40 mechanical forces on bone healing have been previously reviewed 5,6 ; mechanical loading is likely to 41 depend on frequency 7 and cycle number 8 , influences mesenchymal stem cell differentiation 9 , and has a 42 role in guiding healing towards primary and secondary bone healing pathways 10 . In contrast, 43 mechanical loading and relative motion of fragments showed little to no benefit in other studies 11,12 , 44 though the timing of the changes of the mechanical environment in relation to healing phases could 45 also play a role 13 . This further highlights the need to understand the mechanical environment in and 46 around bone defects during healing. This mechanical environment of bone includes the relevant 47 surrounding hard and soft tissues, and their interaction as subject to Newtonian mechanics, which 48 allows computational exploitation for assessing the loading history and relationship to morphological 49 changes 14 . However, the translation of these load-driven bone (re)modelling concepts to highly unique 50 bone defect healing scenarios is lacking. Another current challenge that arises is how to determine the 51 force that needs to be applied in a subject-specific manner in order to have a maximal 52 mechanobiological cue without damaging the bone. 53 54 Finite element (FE) analysis is a well-proven approach to understand micro-scale strains and has been 55 previously successfully used to correlate strain and in vivo bone (re)modelling activities in mice 15-17 . 56 In such workflows, bone mechanoregulation can be studied non-invasively by combining imaging and 57 bridge to adjacent surfaces (Fig. 1b), stabilizing the defect against the applied loading and consequent 136 bending moment induced by the defect asymmetry. bridging. Red/yellow: bone volume formed at week 1, 2, 3, and 4. Blue/purple: bone volume resorbed 142 at week 1, 2, 3, and 4. 143

144
Longitudinal assessment of bone defect healing 145 Two volumes were created to differentiate between the initial empty defect space, and the existing 146 surrounding bone. These were the defect centre (DC) and the defect periphery (DP). (Fig. 2a). Within 147 DC (Fig. 2b), a significant interaction was present for bone volume fraction (BV/TV) between time 148 and loading (F(1,4)=8.90, p<0.005)). At week 2, 3 and 4, loading significantly increased BV/TV over 149 controls (p<0.05). BV/TV increased significantly over time in both loading and control groups 150 (p<0.005). In the DP, both time (F(1,4)=164.7, p<0.005)) and the loading (F(1,4)=5.34, p=0.041) had 151 significant overall effects on BV/TV (Fig. 2c). 152 the rtFE loading from week 2 compared to controls. (c) BV/TV within DP was also found to be 156 influenced, but not to the same magnitude or extent as the DC. 157

158
Overall, loading (F(1,3)=12.6, p<0.001)) and time (F(1,3)=23.2, p<0.001) had a significant effect on 159 the DC bone formation rate (BFR/DC). Loading significantly increased BFR/DC over controls 160 between weeks 1-2 (F(1,3)=5.83, p=0.013), and weeks 3-4 (F(1,3)=4.39, p=0.042), though did not 161 reach significance for weeks 2-3 (F(1,3)=2.05, p=0.159) (Fig. 3a). Loading did not largely affect the 162 DP bone formation rate (BFR/DP) or DP bone resorption rate (BRR/DP) at any time (Fig. 3b). Also, 163 overall, loading did not have a large effect on BRR/DC or BRR/DP compared to control mice (Fig. 3). 164 influenced the BFR/DC compared to controls, and reached significance at postoperative weeks 2 and 167 4, while loading did not appear to influence BRR/DC compared to controls. (b) Loading did not 168 significantly influence either BFR/DP or BRR/DP compared to controls, at any time interval. 169 170 Mechanics of bone volume changes 171 After 4 weeks, none of the vertebrae had recovered their pre-surgery axial stiffness based on the 172 applied forces. From week 2 onwards, the treatment had a significant positive effect on FE-calculated 173 stiffness (p<0.05, Supplementary Fig. S2a). There was also a significant positive Pearson's correlation 174 between BV/TV in the DC, and FE-calculated normalized stiffness (r=0.907, n=51, p<0.005, 175 Supplementary Fig. S2b), while in the DP the correlation to normalized stiffness was significant, but 176 lower (r=0.809, n=51, p<0.005). A pattern was noted in the probability of formation, quiescence, or 177 resorption events within the combined DC and DP regions; they were largely related to effective strain 178 (EFF). The effective strain as a percentile of the max effective strain (EFF/EFFMAX) was used to find 179 the conditional probability of a (re)modelling event occurring due to strain. During the first week of 180 healing, bone formed with a random conditional probability (cp(x) = 33.3%). However, in subsequent 181 weeks, bone formed in the upper half of the strain field (cpF(x > 50 %) > 33.33 %). The probability 182 for resorption was the highest within the first 30% of occurred strain leading to a small strain window 183 where bone was predominantly quiescent between 30% and 50% (Fig. 4a). This pattern also occurred 184 in control animals ( Fig. 4b) where physiological strains were simulated with FE analysis (Fig. 5). 185 When comparing the loading and control cases, the curvature of the formation, resorption and 186 quiescence profiles in the loading cases had steeper and more pronounced curve sections compared to 187 the control cases. There was also a small shift towards lower EFF/EFFMAX being formative with 188 loading.

Figure 4. Conditional probabilities of formation and resorption events within the combined DC 191
and DP regions. Effective strain as a ratio of the maximum effective strain value (EFF/EFFMAX)). (a) 192 Higher ratios of EFF/EFFMAX led to formation activity in treatment groups, and lower ratios led to 193 resorption activity, regardless whether externally loaded in the treatment group, or in (b) control 194 animals with an assumed axial strain. While the influence of mechanical loading on bone (re)modelling is known, implementing this to 202 defect healing creates other challenges. These arise as two aspects, being the changed mechanical 203 environment due to the initial defect, and its unknown healing thereafter. 204

205
This study showed that by introducing an rtFE approach to an existing loading set-up 22 , bone defect 206 healing could be significantly improved over no-treatment controls. Importantly, this approach 207 succeeded in avoiding any incidence of fracture due to overloading, and in principle, homogenized 208 strain across different defect shapes, sizes, and healing progressions. Adaption of the loading within 209 the mechanical environment is not novel in itself, with mixed reports on the effectiveness of 210 dynamization 12,32 , in which a change in stiffness of external fixators is adapted over healing periods. 211 Adding to these existing concepts, this study estimated individualized loads to apply during healing; 212 the rtFE approach allows much greater accuracy in the control of strain, as opposed to generic or pre-213 determined adaptive regimes. 214

215
Loading started at two days after the defect was created, at which time the mice were still relatively 216 young at fourteen weeks old. The bone response to loading has previously been shown to have some 217 age-dependency on mice around this age, where six week old mice had a more exaggerated response 218 to loading compared to ten and sixteen week old mice 33 . Hence, it could be questioned whether the 219 positive effect found in this study would also be repeated in older mouse. In this regard, it has been 220 reported that after sixteen weeks old, aging has less of an influence on the response to loading 30 . At 221 fourteen weeks old, the mice in this study are near the border of this apparent age threshold. As such, 222 prior studies suggest the positive effect of this rtFE loading could be beneficial into adulthood and 223 beyond, though further studies would be needed to confirm this. 224

225
The timing of loading after bone injury is a topic of debate. In general, loading is known to influence 226 cell activity and function due to tissue deformation and fluid flow 34 . Furthermore, it influences both 227 spatial and temporal biological responses at multiple scales 35,36 . This aspect is more important in the 228 case of bone healing, where multiple overlapping phases exist. Early loading from two days onwards 229 has been reported less effective than delayed loading at two weeks 37 . However, it has also been 230 reported that early cyclic loading may increase oxygen transport to the defect region 38 , promoting the 231 longer term regeneration response. As differences in bone parameters were already noted after two 232 weeks in this study, it seems that early loading of three sessions per week is potentially more effective 233 than delayed loading, especially as subject-specific adaptive loading will be in an efficient range 234 without risking damage to early callus structures. 235

236
The new bone that formed within the defect did not appear to be simply recreating the pre-defect bone 237 structure, but to be forming based on other factors. During the defect healing, strain appeared to guide 238 bone (re)modelling activities, more so than an inherent sense of prior bone anatomy. Bone appeared to 239 form in compensation for the asymmetry of the defected bone, without cortical bridging, but with 240 densely arranged trabecular bone within the marrow cavity in a conical, V-like shape. This formation 241 pattern supports the idea of a univerisal cellular mechanobiological response regardless of the location 242 within the body, or the presence of fixation implants. This aligns with concepts of bone resorption 243 being caused by either disuse or stress-shielding, where the cells respond to the loads they experience 244 depending on their mechanical state, and not of the cause of the change. Further, the axial mechanical 245 stiffness over the four weeks, as assessed with FE, did not recover to its pre-defect strength, 246 suggesting that healing would continue into the future. The model was newly developed, and the level 247 of impairment over longer time periods would require further characterization of the model. As this 248 study duration was only four weeks, it is unclear whether (re)modelling would eventually result in a 249 structure similar to the native vertebra over time, or even have healed completely given more time. Regardless, bone formation during early healing was related to and guided by bone strain (Fig. 5). 251 This applied to both the loaded and the control group, which could be exposed to physiological strain 252 that may peak at 4N 14 . Strain-induced bone (re)modelling principles have been previously recognized 15,16 , in which loading favors formation over resorption 17 . This principle is evident 254 elsewhere, where in a study of mouse femoral fractures, fracture site remodeling after three weeks has 255 shown to be consistent with previously considered remodeling theories 39 . However, such comparisons 256 to other locations within the body can be confounded by non-biomechanical factors. For example, the 257 difference in the presence of bone marrow and stem cells, cancellous and cortical bone ratios, and the 258 large muscles surrounding the femur that supply blood and cytokines. The conditional probabilities of 259 formation, quiescence and resorption (Fig. 4)   This study did not compare subject-specific adaptive loading group to any traditional, non-adaptive, 282 non-subject-specific loading group. Mechanical loading is well-established to enhance bone healing, 283 and this has been extensively demonstrated in a variety of animal models using various loading 284 modalities. But many questions remain on how to implement this knowledge into practice, where 285 defects and their healing progression can vary widely. This study developed and implemented an 286 objective 3-dimensional imaging and analysis method to assess a defect and its healing, and 287 demonstrated how this could be linked to a known effective loading regime while reducing secondary 288 fracture risks ( Supplementary Fig. S1). This overall approach attempts to foresee technological 289 progress and tools that could be more reliable than, for example, subjective grading scales of fractures 290 based on 2-dimensional imaging with subsequent loading based on this assigned grading scale. In this 291 regard, this study does not provide evidence that the complex, objective methodology presented in this 292 study provides improved outcomes compared to a simpler subjective analysis and/or non-adaptive 293 loading regime. Future studies could be designed to investigate if such a benefit truly exists between 294 these approaches. In principle, though, a leaning towards objective, (semi-)quantitative analyses have 295 historically prevailed over subjective, qualitative analyses, and this study attempts to follow this path. 296 As discussed, assumptions and simplifications created several limitations to this study which cover 297 both animal and computational aspects. For the animal aspects, the mice were relatively young, the 298 healing was not completed within the four weeks, and three mice were excluded which reduced study 299 power. Defects were created using relatively basic tools, and while this provides simplicity and an 300 ability to apply the model, it also introduces some variability in the defect volume across animals. 301 However, this is factored for in the BV/TV calculation. The pinning of the adjacent vertebra also 302 lacked certain control in position and angle, which creates unknowns in how the defect vertebra is 303 loaded, and the modes it vibrates in considering the semi-constrained nature. Future studies could be 304 extended past four weeks in older mice, improve the defect precision and repeatability, and increase 305 pinning control. As for the computational aspects, one of the greatest limitations is the simplifications 306 in the FE model. Firstly, the materials and models are linear elastic, which would not capture any non-307 linear behaviour of vibration or visco-elastic effects. Material properties of the modelled discs were that of bone, and all surrounding and void voxels not designated as bone were assigned a Young's 309 modulus of 3MPa, including what would be muscle or air. As such, both the disc and non-bone 310 regions are therefore stiffer than reality. Many of these computational limitations relate to the micro-311 FE solver, ParOSol 44 ; however, these compromises enabled the efficient use of running the FE 312 remotely on a supercomputer, which created the possibility for near real-time results. While 313 computational power is itself a limitation, future studies could explore methods to more accurately 314 model such a dynamic system, and further validate the simplifications and assumptions used. 315

316
To translate this approach to patients, not only would further research and development be needed, but 317 that technological advancements to continue. Most obviously, the computational and hardware 318 technology used in this study is not currently available to clinicians. Secondly, any patient-specific 319 solutions within the FE-realm require computational assumptions to be made, which requires further 320 expertise when applying case-by-case. These current limitations in translatability will likely be 321 overcome as technology develops, this study demonstrates the possibilities the future research can 322 strive towards, once technology and methods inevitably catch up for use in the larger-scales needed for 323 humans. 324

325
In conclusion, individualized real-time adaptive loading can be achieved through a combination of 326 micro-CT imaging, followed immediately by FE-solved strain distribution, and finally rescaling and 327 application of a cyclic loading force accordingly. Further investigation is needed to compare this to 328 traditional non-adaptive methods. This rtFE approach is highly relevant for clinical scenarios where 329 bone fractures and their healing progression are unique. This approach optimizes loading intensity, 330 and has the potential to reduce the risk of re-fracture or ineffective mechanical loading, thus improving 331 the healing of bone defects. 332

334
Study design and surgery 335 Approval was obtained for the animal experiments from the cantonal ethics committee from the 336 Kantonales Veterinäramt Zurich (Zurich, Switzerland, ZH029/18) prior to the study, and all 337 experiments were performed in accordance with Swiss animal welfare act and ordinance, and 338 ARRIVE guidelines. The study included two groups: an rtFE loading group, that were adaptively 339 loaded (3.2-5.5N, 10Hz, 5 minutes, 3000 cycles), and; a control group, that received sham loading 340 (0N) and similar handling. Groups were allocated by block randomization within a larger study, with 341 sample sizes estimated from previous similar research within the laboratory 15 , in which 2 groups for a 342 repeated (4) measures ANOVA using G*Power (β=0.8, α=0.05, f =0.7, number measurements=4, 343 correlation=0.8) estimated a total sample size of 16 (n=8 per group). All surgical, scanning and 344 loading procedures were performed under isoflurane anesthesia (induction 5%, maintenance 1-2%, in 345 O2). To be able to apply loading, three weeks prior to defect surgery, stainless steel pins (Fine Science 346 Tools, Heidelberg, Germany) were inserted in the fifth and seventh caudal vertebrae under 347 fluoroscopic control, as previously described (4). Perioperative analgesia (25 mg/L, Tramal®, 348 Gruenenthal GmbH, Aachen, Germany) was delivered via the drinking water for pre-emptive pain 349 relief two days prior to the defect surgery, and for three days post-surgery. All surgeries were The reconstructed images were Gaussian filtered (sigma 1.2, support 1) to reduce noise, and 365 thresholded to assign material properties. Voxels within the threshold range from 395 mgHA/cm 3 to 366 745 mgHA/cm 3 were regarded as bone. This bone was assigned isotropic linear elastic material 367 properties with a Young's modulus between 4 GPa and 12.8 GPa, in threshold steps of 25 mgHA/cm 3 , 368 in proportion with their density 45 . The Young's modulus of soft tissue was set to 3 MPa for values 369 lower than 395 mgHA/cm 3 . The Poisson's ratio was set to 0.3. Vertebra geometry was aligned to the 370 principal axis of the coordinate system in the z-direction. To achieve an even force distribution across 371 the bone and to counter numerical issues with the solver, discs of 1.68 mm diameter with a Young's 372 modulus of 12.8 GPa were added at the distal and proximal ends of the vertebrae, similar to previously 373 used 17 . The outer surface of the distal disc was fixated using Dirichlet boundary conditions. The outer 374 surface of the proximal disc was deformed by 1% by applying a pure compressive force in the normal 375 (z) direction. Each model was solved using the ParOSol solver 44 , running at the Swiss National 376 Supercomputing Centre (CSCS, Lugano, Switzerland) with 64 CPUs, taking less than 2 minutes in 377 computing time. Effective strain was used as the output measurement for its ability to capture 378 inhomogeneity in newly formed bone tissue. The dynamic character of the loading was not considered, 379 in alignment with previous studies reporting that static simulations capture the main features and  bridging. Red/yellow: bone volume formed at week 1, 2, 3, and 4. Blue/purple: bone volume resorbed 614 at week 1, 2, 3, and 4. 615 the rtFE loading from week 2 compared to controls. (c) BV/TV within DP was also found to be 618 influenced, but not to the same magnitude or extent as the DC. 619 influenced the BFR/DC compared to controls, and reached significance at postoperative weeks 2 and 621 4, while loading did not appear to influence BRR/DC compared to controls. (b) Loading did not 622 significantly influence either BFR/DP or BRR/DP compared to controls, at any time interval. 623 and DP regions. Effective strain as a ratio of the maximum effective strain value (EFF/EFFMAX)). (a) 625 Higher ratios of EFF/EFFMAX led to formation activity in treatment groups, and lower ratios led to 626 resorption activity, regardless whether externally loaded in the treatment group, or in (b) control 627 animals with an assumed axial strain. 628 Figure 5. Representative sample depicting formation and resorption relationship with effective 629 strain. Regions of higher effective strain tended to formation, while lower effective strain tended to 630 resorption. 631