Mechano-regulation of bone adaptation is controlled by the local in vivo environment and logarithmically dependent on loading frequency

It is well established that cyclic, but not static, mechanical loading has anabolic effects on bone. However, the function describing the relationship between the loading frequency and the amount of bone adaptation remains unclear. Using a combined experimental and computational approach, this study aimed to investigate whether bone mechano-regulation is controlled by mechanical signals in the local in vivo environment and dependent on loading frequency. Specifically, by combining in vivo micro-computed tomography (micro-CT) imaging with micro-finite element (micro-FE) analysis, we monitored the changes in microstructural as well as the mechanical in vivo environment (strain energy density (SED) and SED gradient) of mouse caudal vertebrae over 4 weeks of either cyclic loading at varying frequencies of 2Hz, 5Hz, or 10Hz, respectively or static loading. Higher values of SED and SED gradient on the local tissue level led to an increased probability of bone formation and a decreased probability of bone resorption. In all loading groups, the SED gradient was superior in the determination of local bone formation and resorption events as compared to SED. Cyclic loading induced positive net remodeling rates when compared to sham and static loading, mainly due to an increase in mineralizing surface and a decrease in eroded surface. Consequently, bone volume fraction increased over time in 2Hz, 5Hz and 10Hz (+15%, +21% and +24%, p<0.0001), while static loading led to a decrease in bone volume fraction (−9%, p≤0.001). Furthermore, regression analysis revealed a logarithmic relationship between loading frequency and the net change in bone volume fraction over the four week observation period (R2=0.74). In conclusion, these results suggest that bone adaptation is regulated by mechanical signals in the local in vivo environment and furthermore, that mechano-regulation is logarithmically dependent on loading frequency with frequencies below a certain threshold having catabolic effects, and those above anabolic effects. This study thereby provides valuable insights towards a better understanding of the mechanical signals influencing bone formation and resorption in the local in vivo environment.

as the mechanical in vivo environment (strain energy density (SED) and SED gradient) of 23 mouse caudal vertebrae over 4 weeks of either cyclic loading at varying frequencies of 2Hz, 24 5Hz, or 10Hz, respectively or static loading. Higher values of SED and SED gradient on the 25 local tissue level led to an increased probability of bone formation and a decreased probability 26 of bone resorption. In all loading groups, the SED gradient was superior in the determination 27 of local bone formation and resorption events as compared to SED. Cyclic loading induced 28 positive net remodeling rates when compared to sham and static loading, mainly due to an 29 increase in mineralizing surface and a decrease in eroded surface. Consequently, bone volume 30 fraction increased over time in 2Hz, 5Hz and 10Hz (+15%, +21% and +24%, p<0.0001), while 31 static loading led to a decrease in bone volume fraction (-9%, p≤0.001). Furthermore, regression 32 analysis revealed a logarithmic relationship between loading frequency and the net change in 33 bone volume fraction over the four week observation period (R 2 =0.74). In conclusion, these 34 results suggest that bone adaptation is regulated by mechanical signals in the local in vivo 35 environment and furthermore, that mechano-regulation is logarithmically dependent on loading 36 frequency with frequencies below a certain threshold having catabolic effects, and those above 37 anabolic effects. This study thereby provides valuable insights towards a better understanding 38 Introduction 43 It is well established that cyclic, but not static loading has anabolic effects on bone [1][2][3][4]. This 44 clear-cut discrepancy in osteogenic responses to both loading patterns highlights the key role 45 of loading frequency in mechano-regulation of bone remodeling -the coordinated process by 46 which bone is continuously formed and resorbed. Yet, the exact relationship between loading 47 123 Figure 1 shows the relative changes in trabecular bone morphometric parameters over the 4-124 week loading period for the different loading groups. BV/TV developed differently over time 125 between the loading groups (interaction effect, p<0.0001). Compared to the sham-loaded group, 126 which showed no change in BV/TV over time (-6%, p>0.05), cyclic loading at all frequencies 127 (2Hz, 5Hz and 10Hz) led to a dose-response increase in BV/TV with higher frequencies 128 resulting in higher increases in BV/TV ( Fig 1A). Herein, the 5Hz and 10Hz groups showed a 129 significant increase compared to baseline already 2 weeks after the start of loading (p≤0.001 130 and p<0.0001), while the 2Hz group showed a significant increase relative to baseline only after 131 three weeks (p≤0.001). At the end of the 4-week loading regime, these groups showed a 15%, 132 21% and 24% higher BV/TV relative to baseline (p<0.0001 for 2Hz, 5Hz and 10Hz). Static 133 loading on the other hand, had catabolic effects resulting in significantly lower BV/TV (-9%, 134 p≤0.01) at the last time point relative to baseline. In line with the changes in BV/TV, Tb.Th 135 developed differently over time between the loading groups (interaction effect, p<0.0001, Fig  136   1B). By the end of the 4-week loading intervention, all cyclic loading groups showed significant 137 increases in Tb.Th (p<0.0001), which was not observed in the static and sham-loaded groups 138 (p>0.05). Although the number of trabeculae (Tb.N) decreased and trabecular separation 139 increased (Tb.Sp) over time (Fig 1C,D,p<0.001), no relative differences were observed 140 between the groups (p>0.05). These results thus suggest that increases in BV/TV due to cyclic 141 loading were mainly driven by thickening of the trabeculae rather than by the inhibition of the 142 reduction in the number of trabeculae. 143 By plotting the relative changes in BV/TV as a function of loading frequency, regression 144 analysis revealed a logarithmic relationship between bone adaptation and loading frequency as 145 a best fit to the data (R 2 =0.74, Fig 1F) with loading frequencies above 0.36Hz±0.08 having 146 anabolic effects, and frequencies below this threshold having catabolic effects. Although there 147 were no significant differences between the cyclic loading groups, loading at 10Hz had the 148 earliest and largest anabolic effects compared to the other frequencies. value for main effect of group determined by one-way ANOVA, **** p<0.0001 denotes 159 significant difference between groups determined by post hoc Tukey's multiple comparisons 160 test). 161 Aside from providing information on changes in bone structural parameters over time, in vivo 162 micro-CT also provided the possibility to assess dynamic bone formation and resorption 163 activities such as bone formation/resorption rate (BFR/BRR), mineral apposition/resorption 164 rate (MAR/MRR) and mineralizing/eroded surface (MS/ES) [18]. The net remodeling rate 165 (BFR-BRR), which gives an indication whether there was overall bone gain (i.e., BFR-BRR>0) 166 or loss (i.e., BFR-BRR<0) occurring within the trabecular compartment, tended to develop 167 differently between groups (p=0.056). Compared to the static and sham-loaded groups, which 168 had an overall negative remodeling balance, the 2Hz, 5Hz and 10Hz had an overall positive 169 remodeling balance (p≤0.01, p≤0.001 and p<0.0001, Fig 2A). The net remodeling rate did not 170 significantly change over time. When bone formation and resorption rates were analyzed 171 separately, the main differences in the cyclic loading groups were in the reduced BRR as 172 compared to the sham and static groups. While BFR did not significantly differ between groups 173 (p>0.05, Fig 2B), BRR was 35% (p<0.01), 50% (p<0.0001) and 44% (p<0.0001) lower in the 174 2Hz, 5Hz, 10Hz groups, respectively, compared to the sham-loaded group ( Fig 2C). The static 175 group on the other hand had a similar BRR (-2%, p>0.05) as the sham-loaded group ( Fig 2C). 176 A difference between the cyclic and static loading groups was also apparent when investigating 177 the surfaces of formation (mineralized surface, MS, interaction effect p<0.05) and resorption 178 (eroded surface, ES, interaction effect p≥0.05) sites with the cyclic loading groups having a 179 higher MS and lower ES compared to the static and sham-loaded groups ( Fig 2D). On average, 180 formation sites occupied 2, 2.5 and 2.6 more surfaces than resorption sites for the 2Hz, 5Hz and 181 10Hz groups, and only 1.4 times more for the control and static groups, respectively. Furthermore, the 2Hz, 5Hz and 10Hz groups had a 18% (p=0.0078), 25% (p=0.0007) and 26% 197 (p>0.0001) higher mineralized surface (MS) and a 22% (p<0.0001), 32% (p<0.0001) and 26% 198 (p<0.0001) lower eroded surface (ES) compared to the sham-loaded group, while the static 199 group had similar MS and ES compared to sham-loading (p>0.05, Fig 2E-F). The mineral 200 apposition and resorption rates (MAR and MRR), which represent the thicknesses of formation 201 and resorption packages, respectively, did not develop differently between groups (interaction 202 effects p=0.586 and p=0.459). Furthermore, the MAR and MRR were similar between groups 203 (p>0.05), thus suggesting that they are not affected by loading (Fig 2G-I). This indicates that 204 cyclic loading had a greater effect on surface than on thickness of formation as well as 205 resorption sites. 206

Bone adaptation to load is controlled by mechanical signals in the local in vivo 207
environment 208 In order to assess whether bone remodeling events -namely formation, quiescence (i.e., where 209 no remodeling occurred) and resorption -can be linked to the corresponding mechanical signals 210 in the local in vivo environment, we performed micro-finite element (micro-FE) analysis to 211 calculate the strain distribution within the tissue. As deformation (direct cell strain) and well as the corresponding maps of SED ( Fig 3B) and ∇SED (Fig 3C). From this qualitative 218 analysis, it is apparent that bone resorption occurs at sites of lower SED and ∇SED, respectively, 219 whereas bone formation occurs at sites of higher SED and ∇SED (Fig 3). images showing sites of bone formation (orange), quiescence (grey) and resorption (purple). 225 Corresponding map of the (B) strain energy density (SED) and (C) gradient thereof (∇SED) 226 showing sites of higher (red) and lower (blue) SED/∇SED values obtained by micro-finite 227 element (micro-FE) analysis. 228 To establish a quantitative description of the mechano-regulation of bone remodeling, we 229 calculated the conditional probabilities for a given remodeling event to occur as a function of 230 the mechanical stimuli, also known as remodeling rules [15]. Figure 4 shows the conditional 231 probability curves for formation (orange), quiescence (grey) or resorption (purple) to occur at 232 a given value of SED (Fig 4A,C,E) or ∇SED (Fig 4B,D,F) for the different groups averaged 233 over all time points. For all groups, the conditional probability for bone formation to occur was 234 higher at higher values of SED and ∇SED, respectively (SED/SEDmax > 0.18) whereas bone 235 resorption was more likely to occur at lower values (SED/SEDmax < 0.18). The probability 236 curves for all groups were fit by exponential functions (Table S1), of which the coefficients 237 provide information on the functioning of the mechanosensory system as described previously 238 [15]. When comparing the slopes of the formation probability curves (parameter a, Fig 4A,B  239 and Table S1), which can be interpreted as the mechanical sensitivity of the system, there was 240 a gradual increase of the mechanical sensitivity with increasing frequency with the 10Hz group 241 showing the highest mechanical sensitivity (aSED = 0.217, aSEDgrad = 0.316). For the resorption 242 probability curves (Fig 4E,F and Table S1), the 5Hz and 10Hz groups showed similar 243 mechanical sensitivity to SED (aSED = 0.284), while the 5Hz group showed highest sensitivity 244 to ∇SED (aSEDgrad = 0.264 compared to aSEDgrad = 0.252 in 10Hz group). The probability of the 245 quiescence however, was not influenced by loading frequency (Fig 4C,D). When comparing 246 between SED and ∇SED as mechanical stimuli driving bone remodeling events, it seems that in 247 all groups, formation was more sensitive to ∇SED shown by the higher slopes (aSED < agradSED) 248 of the probability curves (Fig 4A,B and Table S1). In contrast, resorption seemed to be more To better compare the modeling performance of SED versus ∇SED for the prediction of bone 257 remodeling events, an area under the receiver operator characteristic curve (AUC) approach 258 was used (Fig 5). For all groups, the AUC values for formation (for all groups p<0.0001, Fig  259   5A) and resorption (for all groups p<0.05 except for 5Hz p<0.10, Fig 6C) events were higher 260 for the ∇SED compared to SED. No difference between SED and ∇SED was observed for 261 quiescence (Fig 5B). These results suggest that ∇SED has a better modeling performance 262 compared to SED for determining the probability of bone formation and resorption events. While static loading had catabolic effects, cyclic loading at 2Hz, 5Hz and 10Hz had anabolic 277 effects on trabecular bone. In line with previous studies using the tail loading model [17,26], 278 cyclic loading over four weeks led to an increase in BV/TV, which was driven by the thickening 279 of individual trabeculae rather than a prevention of loss in trabecular number. Furthermore, by 280 registering consecutive time-lapsed images onto one-another, we were able to quantify both 281 bone formation as well as bone resorption activities in three dimensions [26], which to the best 282 of our knowledge, has not yet been used to assess the effects of static loading regimes. 283 Specifically, we showed that cyclic loading mainly affects the surfaces of the bone formation 284 and resorption sites (MS and ES), rather than the thickness of these remodeling packets (MAR 285 and MRR). In agreement with previous studies [18,26], these results suggest that cyclic loading 286 promotes osteoblast recruitment, while simultaneously inhibiting osteoclast recruitment. 287 Ultimately, cyclic loading results in larger mineralized surfaces and smaller eroded surfaces 288 while keeping the thickness of the remodeling packets constant. 289 Notably, this study showed a logarithmic relationship between loading frequency and load-290 induced bone adaptation with frequencies above a certain threshold having anabolic effects and 291 those below having catabolic effects. That cyclic, but not static loading, has anabolic effects on 292 cortical bone has been shown in various animal models including rabbits [2], turkeys [1] and 293 rats [3,4]. However, to the best of our knowledge, the effect of static loading has not yet been 294 assessed in trabecular bone in mice. In line with the existence of a frequency threshold 295 (0.36Hz±0.08) to elicit anabolic responses as demonstrated in this study, Turner et al. found 296 that bone formation rate in rat tibiae only increased with frequencies above 0.5Hz, followed by 297 a dose-response increase up to 2Hz [5]. Using a similar design as our study, Warden et al. 298 showed increased cortical bone adaptation with increasing loading frequencies up to 5 to 10Hz 299 with no additional benefits beyond 10Hz [10]. In a theoretical study, Kameo et al. furthermore 300 showed similar results by subjecting individual trabeculae to uniaxial loading at frequencies 301 ranging from 1 to 20Hz [12]. Although one would expect higher loading frequencies to lead to 302 higher cellular stimulation and a consequent greater anabolic response, it has been suggested 303 that frequencies above a certain threshold (10Hz) reduce the efficiency of fluid flow through 304 the LCN, thus resulting in inefficient mechanotransduction [10,29]. More recently, by 305 monitoring Ca 2+ signaling in living animals, Lewis et al. have shown that osteocyte recruitment 306 was strongly influenced by loading frequency [30]. Another physiological system, for which 307 the relationship between frequency and mechanotransduction is widely studied, is the inner ear 308 [31,32]. Hair cells, the cells responsible for transducing mechanical forces originating from 309 acoustic waves to neural signals, are sensitive to frequency [31,33]. Furthermore, the sensitivity 310 of the ear varies with the frequency of sound waves resulting in a limited range of frequencies 311 that can be perceived. Hence, drawing an analogy to the theory of sound pressure level, which 312 also displays logarithmic laws [34], it is possible that bone's response to frequency is similar 313 to the perception of sound in human hearing. 314 One limitation of this study was that loading at low (1Hz) and higher (>10Hz) frequencies was 315 not assessed. Furthermore, as the strain magnitude and duration of individual loading bouts 316 were the same for all loading groups, the number of cycles and strain rate differed between the 317 different loading groups. From this study design, it therefore remains impossible to know 318 whether the number of cycles or the loading frequency are the main factors driving load-induced 319 bone adaptation. Hence, whether bone's osteogenic response to loading is indeed limited to a 320 specific range of frequencies, below and above which bone becomes less osteogenic, requires 321 further in vivo experiments. 322 Using the combined approach of time-lapsed in vivo micro-CT imaging and micro-FE analysis, 323 we showed that bone remodeling activities were correlated to the local mechanical environment 324 at the tissue level. In agreement with previous studies [15,17], bone formation was more likely 325 to occur at sites of higher SED whereas bone resorption was more likely to occur at sites of 326 lower SED. Furthermore, compared to static loading, cyclic loading decreased the probability 327 of non-targeted bone remodeling, which led to an increase in bone formation and a decrease in 328 bone resorption. In addition, we showed that the SED gradient was better at predicting bone 20 formation and resorption events compared to SED. That the SED gradient, a measure of fluid 330 flow through the LCN, can improve predictions of remodeling events compared to SED, a 331 measure of direct cell strain, has been suggested previously [16]. Furthermore, as the SED 332 gradient encompasses the neighboring SED voxels, it provides information of a broader 333 mechanical environment, which could explain the higher modeling performance observed with 334 the SED gradient compared to SED. An additional limitation of this study was that the micro-335 FE analysis did not take into account the component of frequency. Although our approach 336 enabled us to link bone remodeling events to mechanical environments in vivo at the local level, 337 the addition of theoretical models that incorporate cellular mechanosensing and intercellular 338 communication [12,[35][36][37] will be highly useful to improve our understanding of the 339 relationship between loading frequency and bone adaptation across multiple scales. 340 In conclusion, these results suggest that bone adaptation is regulated by mechanical signals in 341 the local in vivo environment and furthermore, that mechano-regulation is logarithmically 342 dependent on loading frequency with frequencies below a certain threshold having catabolic 343 effects, and those above anabolic effects. This study thereby provides valuable insights towards 344 a better understanding of the mechanical signals influencing bone formation and resorption in 345 the local in vivo environment. 346 347

Study Design 349
To investigate the effect of loading frequency on mouse caudal vertebrae, 11-week old female 350 C57BL/6J mice were purchased (Charles River Laboratories, France) and housed at the ETH 351 Phenomics Center (12h:12h light-dark cycle, maintenance feed and water ad libitum, three to 352 five animals/cage) for one week. To enable mechanical loading of the 6 th caudal vertebrae 353 (CV6), stainless steel pins (Fine Science Tools, Heidelberg, Germany) were inserted into the 354 fifth and seventh caudal vertebrae of all mice at 12 weeks of age. After three weeks of recovery, 355 the mice received either sham (0N), 8N static or 8N cyclic loading with frequencies of 2Hz, 356 5Hz, or 10Hz and were scanned weekly using in vivo micro-CT. All procedures were performed 357 under isoflurane anaesthesia (induction/maintenance: 5%/1-2% isoflurane/oxygen). All mouse 358 experiments described in the present study were carried out in strict accordance with the 359 recommendations and regulations in the Animal Welfare Ordinance (TSchV 455.1) of the Swiss 360 Federal Food Safety and Veterinary Office (license number 262/2016). 361

Mechanical loading 362
The loading regime was performed for five minutes, three times per week over 4 weeks as 363 described previously [14]. For the cyclic loading groups, sinusoidally varying forces (8N 364 amplitude) were applied at 2Hz, 5Hz or 10Hz resulting in cycle numbers of 600, 1500 and 3000, 365 respectively. For the static loading group, the force was maintained at 8N during the five 366 minutes. For the sham-loaded group, the tails were fixed in the loading device for five minutes, 367 but no loading was applied (0N). 368

Micro-CT imaging and analysis 369
In vivo micro-CT (vivaCT 40, Scanco Medical AG, isotropic nominal resolution: 10.5 µm; 55 370 kVp, 145 µA, 350 ms integration time, 500 projections per 180°, scan duration ca. 15 min, radiation dose per scan ca. 640 mGy) images of the CV6 were acquired every week. Micro-CT 372 data was processed and standard bone microstructural parameters were calculated in trabecular, 373 cortical and whole bone by using automatically selected masks for these regions as described 374 previously [26]. To calculate dynamic morphometric parameters, micro-CT images from 375 consecutive time-points were registered onto one another. The voxels present only at the initial 376 time point were considered resorbed whereas voxels present only at the later time point were 377 considered formed. Voxels that were present at both time points were considered as quiescent 378 bone. By overlaying the images, morphometrical analysis of bone formation and resorption 379 sites within the trabecular region allowed calculations of bone formation rate (BFR), bone 380 resorption rate (BRR), mineral apposition rate (MAR), mineral resorption rate (MRR), 381 mineralizing surface (MS) and eroded surface (ES) [18]. 382

Micro-finite element (micro-FE) analysis 383
For each mouse at each time point, segmented image data was converted to 3D micro-FE 384 models, with additional voxels added to the proximal and distal ends of the vertebrae mimicking 385 intervertebral discs. All voxels were converted to 8 node hexahedral elements and assigned a 386 Young's modulus of 14.8 GPa and a Poisson's ratio of 0.3 [14]. The bone was assumed to have 387 linear elastic behaviour, which allowed for static loading in the micro-FE analysis [38]. The top 388 was displaced by 1% of the length in z-direction (longitudinal axis), while the bottom was 389 constrained in all directions. The micro-FE model was solved using a micro-FE solver 390 (ParOSol). The results were then rescaled to an applied force of 8N for the loaded groups and 391 4N (physiological loading) for the sham-loaded group (0N) as described previously [39]. 392

Mechanical environment 393
The mechanical stimuli, which are hypothesized to drive load induced bone adaptation are 394 deformation (direct cell strain) and interstitial fluid flow (shear stress) [27]. As a measure of 395 the mechanical deformation, strain energy density (SED) magnitudes, defined as the increase 396 in energy associated with the tissue deformation per unit volume, were analysed on the bone 397 surface on the marrow-bone interface. Furthermore, based on the assumption that spatial 398 differences in tissue deformation induce fluid flow, the spatial gradient of the SED was 399 analyzed on the marrow side of the marrow-bone interface [28]. The spatial gradients in x, y 400 and z-direction were calculated as follows: 401 , , Where is the SED of a voxel at x, y, z-position i, Nx,y,z the number of voxels in the x,y,z-403 direction and a the nominal resolution. The norm of the gradient vector (∇SED) was used as a 404 quantity for the fluid flow as described previously [16]. 405

∇SED = √ ( ) + ( ) + ( ) 406
The conditional probabilities for a certain remodeling event (formation, quiescence, resorption) 407 to occur at a given value of SED and ∇SED were calculated as described previously [15]. 408 Briefly, the surface SED and ∇SED values were normalized within each animal and 409 measurement by the maximal SED or ∇SED, respectively (chosen as the 99 th percentile of the 410 values present at the surface and in the volume of interest (VOI) in order to remove the variance 411 due to temporal bone adaptation, applied force in FE analysis and individual animals. For each 412 region (formation, quiescence and resorption), a frequency density histogram with 50 bins and 413 equal bin width was created. In order to rule out the dependence on the imbalance between bone 414 formation and resorption, all remodeling events were assumed to have the same occurrence 415 probability (i.e., formation, resorption and quiescent regions were rescaled to have the same 416 amount of voxels). The remodeling probabilities were fitted by exponential functions using 417 non-linear regression analysis. 418 To quantify the modeling performance of SED and ∇SED, respectively, the area under the curve 419 (AUC) of a receiver operating characteristic (ROC) curve was used. The AUC can be defined 420 as the probability that a randomly selected case ("true") will have a higher test result than a 421 randomly selected control ("false") [40]. The ROC curve is a binary classifier, therefore the 422 three different surface regions were analysed separately and only voxels and mechanical 423 quantity values on the bone or marrow surface were used for the classification.