Bone remodelling in the mouse tibia is spatio-temporally modulated by oestrogen deﬁciency and external mechanical loading: A combined in vivo / in silico study

Osteoporosis disrupts the healthy remodelling process in bone and affects its mechanical properties. Me- chanical loading has been shown to be effective in stimulating bone formation to mitigate initial bone loss. However, no study has investigated the effects of repeated mechanical loading, with a pause of one week in between, in the mouse tibia with oestrogen deﬁciency. This study uses a combined experimental and computational approach, through longitudinal monitoring with micro-computed tomography, to evaluate the effects of loading on bone adaptation in the tibiae of ovariectomised (OVX) C57BL/6 mice from 14 to 22 weeks of age. Micro-FE models coupled with bone adaptation algorithms were used to estimate changes in local tissue strains due to OVX and mechanical loading, and to quantify the relationship between local strain and remodelling. The ﬁrst in vivo mechanical loading increased apposition, by 50–150%, while resorption decreased by 50–60%. Both endosteal and periosteal resorption increased despite the second mechanical loading, and periosteal resorption was up to 70% higher than that after the ﬁrst loading. This was found to correlate with an initial decrease in average strain energy density after the ﬁrst loading, which was lower and more localised after the second loading. Predictions of bone adaptation showed that between 50 and 90% of the load-induced bone apposition is linearly strain driven at the organ-level, but resorption is more biologically driven at the local level. The results imply that a systematic increase in peak load or loading rate may be required to achieve a similar bone adaptation rate in speciﬁc regions of interests.


a b s t r a c t
Osteoporosis disrupts the healthy remodelling process in bone and affects its mechanical properties. Mechanical loading has been shown to be effective in stimulating bone formation to mitigate initial bone loss. However, no study has investigated the effects of repeated mechanical loading, with a pause of one week in between, in the mouse tibia with oestrogen deficiency. This study uses a combined experimental and computational approach, through longitudinal monitoring with micro-computed tomography, to evaluate the effects of loading on bone adaptation in the tibiae of ovariectomised (OVX) C57BL/6 mice from 14 to 22 weeks of age. Micro-FE models coupled with bone adaptation algorithms were used to estimate changes in local tissue strains due to OVX and mechanical loading, and to quantify the relationship between local strain and remodelling. The first in vivo mechanical loading increased apposition, by 50-150%, while resorption decreased by 50-60%. Both endosteal and periosteal resorption increased despite the second mechanical loading, and periosteal resorption was up to 70% higher than that after the first loading. This was found to correlate with an initial decrease in average strain energy density after the first loading, which was lower and more localised after the second loading. Predictions of bone adaptation showed that between 50 and 90% of the load-induced bone apposition is linearly strain driven at the organ-level, but resorption is more biologically driven at the local level. The results imply that a systematic increase in peak load or loading rate may be required to achieve a similar bone adaptation rate in specific regions of interests.

Statement of Signiicance
Bone is a dynamic tissue that adapts its architecture and material properties to changes in external loading. Understanding how bone apposition and resorption are affected by several repeated periods of mechanical stimuli is important for the clinical management of osteoporosis. The loadinduced bone changes after the first and second week of mechanical loading are quantified for the first time in the OVX mouse tibia. Using finite element analysis, we also evaluated the extent bone adaptation is strain driven for each period of loading, and whether this occurs at the organ or local level. This combined experimental-computational approach

Introduction
Bone is a dynamic tissue that adapts its structural properties to changes in mechanical demands, by altering its (re)modelling activities [1] . Bone diseases such as osteoporosis and osteoarthritis affect the bone (re)modelling and the bone mechanical properties. In particular, osteoporotic patients suffer from a systemic reduction in bone mineral density (BMD), bone quality and thickness, caused by the disruption of normal bone adaptation mechanisms. Understanding the effects of external loads on bone (re)modelling is fundamental to advance the development of rehabilitation schemes and interventions to treat osteoporosis and other musculoskeletal diseases [ 2 -4 ], in order to lower the socio-economic burden on society.
In vivo loading of the mouse tibia is commonly used for studying the response of bone adaptation to external loading [ 1 , 5 -7 ]. In this animal model, the murine tibia is subjected to dynamic cyclic compressive loads, and micro-computed tomography (micro-CT) scans can be acquired to quantify the changes in structural properties with time. To ascertain the strain distribution in the tibia, strain gauges [ 1 , 5 , 8 ], digital image correlation (DIC) [9] , digital volume correlation (DVC) [10] and/or finite element (FE) models [ 7-9 , 11 , 12 ] are used. Several studies have shown that under the application of the same load, the strain magnitudes induced in the cortical bone were higher in tibiae from old mice than in tibiae from young mice, due to the reduction in cortical thickness and moment of inertia with age [ 5 , 7 ]. Ovariectomy is a model of oestrogen deficiency that induces accelerated bone resorption [13] , resulting in lower bone mass and strength [14] , and is thus used extensively in murine studies for the optimisation of treatments against skeletal diseases [15] . In ovariectomised (OVX) mice, mechanical loading of the caudal vertebra and tibia has reported improvements in trabecular bone volume fraction (BV/TV) and trabecular thickness [ 3 , 16 ]. In the murine tibia, persistent increases in cortical thickness, cortical area, and bone mineral content (BMC) were also found with age [16] . However, OVX reduces the availability of oestrogen receptors to process strain signals [17] , resulting in the additive effect of mechanical loading to be mitigated with time, and was lost by 26 weeks, in murine tibiae that were axially loaded for 4 weeks prior to OVX [18] .
The effect of strain-matched loading on cortical bone formation with age has yielded a plethora of results, with studies reporting a reduction in cortical bone formation between 10 and 78 weeks old mice under 1200 microstrain [19] , or an increase in cortical bone formation between 8 and 20 weeks old mice under 20 0 0 microstrain [1] . Moreover, regional differences in cortical bone adaptation have been reported. Thus micro-FE models built from micro-CT images are useful in understanding how the strain distribution due to the morphology of the tibia might have led to the differences in bone adaptation, and to understand the local changes that ensued from the treatment [ 6 , 8 , 9 , 12 ]. In OVX mice, this combined micro-CT and micro-FE approach has primarily been applied to the caudal vertebra [ 3 , 4 ]. Conducting such studies on the mouse tibia offers several advantages as it is a non-invasive approach to study bone remodelling in a load bearing site, and it could be used to evaluate the effects of physiological and para-physiological loading on bone adaptation [ 12 , 20 ]. Moreover, it could be integrated to higher dimensional scales with gait analysis in mice [21] and to test the effects of biomaterials [22] .
Prediction of bone adaptation has been used in FE models to evaluate the effects of osteoporosis, and/or to understand the contribution of mechanical stimuli on bone (re)modelling [ 3 , 4 , 6 , 9 , 11 ]. This combination of structural FE models and a mechanoregulation algorithm that drives bone changes based on the local mechanical stimuli predicted the changes in BV/TV and trabecular thickness in the murine caudal vertebra, with maximum errors of 4.5% and 8.8%, respectively [ 3 , 23 ]. Moreover, these algorithms were able to predict the local remodelling sites due to mechanical loading in healthy bones with reasonable accuracy, achieving a spatial match of 55% in the caudal vertebrae [23] and a Kendall's trank coefficient of 0.51 in the murine tibiae [6] . The spatial match of remodelling sites following OVX was 47.6% in the caudal vertebrae, but no previous study has investigated the accuracy of using the mechanoregulation algorithm to predict the effects of osteoporosis and subsequent mechanical loading on the properties of the murine tibiae, to determine if the spatial changes are linearly driven. Moreover, in the caudal vertebra study, the same set of bone remodelling parameters was used to predict bone adaptation throughout the period of mechanical loading intervention [ 3 , 4 ]. Their results showed this method performed well only when the rates of change in densitometric properties remained similar after each set of mechanical loading [3] . However, the structural properties of murine tibia have been reported to be affected by ageing and thus may not exhibit a similar response after each week of loading [7] .
A recently developed algorithm has been able to predict cortical bone changes in the murine tibiae under physiological loading, achieving a similar spatial agreement of 49% and 59% in resorption and apposition, respectively [11] . The optimisation of the parameters in this model was achieved by matching the volumetric second moment [11] between each scan instead of matching the average BV/TV in the caudal vertebra mentioned previously [ 3 , 4 ], as the second moments are sensitive in locating regions of bone adaptation [ 7 , 9 ].
The aim of this study was to use a combined experimental and computational approach to evaluate the effects of OVX and two applications of mechanical loading (3 sessions per week, separated by a week in between) on the bone adaptation and on the changes in micromechanical properties in the murine tibiae. Densitometric changes throughout the tibial length and the number of local sites of bone adaptation were quantified between each micro-CT scan. The effects of OVX and loading on the mechanical properties of the mouse tibia were simulated using micro-FE models and a bone remodelling framework [11] , to assess the extent the adaptive response observed with in vivo micro-CT imaging was linearly related to the local stimulus under physiological loading. Physiological loading, rather than the actual mechanical loads [ 5-7 , 9 , 12 , 20 ], were applied in the micro-FE models to reflect the effects of these adaptive response on normal living conditions for the first time.

Experimental in vivo data
Six female C57BL/6 J mice were ovariectomised (OVX) at week 14 of age to induce oestrogen deficiency, a common model for osteoporosis. The mice were scanned every two weeks from week 14 until week 22 using an in vivo micro computed tomography (micro-CT; vivaCT 80, Scanco Medical, Brüttisellen, Switzerland) protocol that allows for the scanning of the whole tibia with minimal effect of radiation on bone remodelling [ 24 , 25 ]. The X-ray source was set at 55 keV and 145 µA, with an isotropic voxel size of 10.4 µm, field of view of 32 mm, integration time of 10 0 ms and 1500/750 samples/projections [ 24 , 25 ]. A third-order polynomial beam hardening correction algorithm based on a 1200 mg HA/cm 3 wedge phantom was applied to all the scans [ 24 , 25 ]. The right tibia of each mouse was fixed between two soft cups and mechanically loaded using an in vivo tibial loading rig and protocol adapted from de Souza et al. [1] , which has been shown to induce lamellar bone adaptation in both cortical and trabecular bone [ 1 , 9 ], and maximum cortical bone formation at 12 N without inducing microdamage in the bone. A 2-12 N trapezoidal load was applied at weeks 19 and 21, 40 cycles a day and 3 times a week on alternate days (ElectroForce BioDynamic 510 0, TA instruments) ( Fig. 1 ). Peak loads of 12 N were held for 0.2 s, and a 10 s interval separated each load cycle. Loading was achieved by superimposing a dynamic load (16 kN/s) on a static 2.0 N preload. This approach was found to induce compressive strains of approximately −150 0 µε and tensile strains of 20 0 0 µε in the mid-shaft of healthy 12 weeks old tibiae [1] . Moreover, a recent study using DVC confirmed that the loading approach would induce higher strains in the proximal metaphysis compared to the mid-shaft [10] . Loading took place on alternate weeks from the micro-CT scanning Please cite this article as: V.S. Cheong, B.C. Roberts   to avoid two anaesthesia sessions in the same week for the welfare of the animals. The experiments were matched in peak force rather than in strain to highlight the effect of the microstructural changes on the local mechanical properties of the mouse tibia after the two separate weeks of loading, and to improve the translatability of the findings in the prescription of loads for rehabilitation exercises. The mice were weighed weekly, before and after OVX surgery, and were 21.2 ± 1.8 g throughout the experimental study.
The study design also consisted of two longitudinal control groups of OVX and wild type (WT) mice (positive and negative controls, respectively), to compare the effects of OVX between weeks 14 and 18, and the effects of loading between weeks 18 and 22; the tibiae were imaged with the same scanning protocol as in this study [13] . A SHAM group was not conducted in line with two of the 3Rs principles to reduce and refine the number of animals used in experiments [26] . Weight gain due to OVX was not evaluated in this study, but in the previous study [13] , where the weight of the mice were 22.6 ± 3.1 g and 20.1 ± 1.8 g in the positive and negative controls, respectively. All experimental procedures were conducted in compliance with the ARRIVE (Animal Research: Reporting of in vivo Experiments) guidelines, the UK Animals (Scientific Procedures) Act 1986 and were approved by the local Research Ethics Committee of the University of Sheffield.

Bone geometry
To compare the changes in bone geometries and to ensure a similar alignment for all images at every time point, the micro-CT images were rigidly registered to a reference bone, as detailed in Lu et al. [ 26 , 27 ]. This registration procedure has been shown to achieve less than 3.5% error, and an intraclass correlation coefficient above 0.8 in local bone mineral content [ 26 , 27 ]. The presence of the fibula affected the rigid registration process and thus the fibula was virtually removed from the images prior to the registration of the bones images (Amira 6.3.0, FEI Visualization Sciences Group, France [Mutual information metric, Lanczos interpolation]). The fibula was removed using a semi-automatic script developed in-house that creates two disconnections near the tibio-Please cite this article as: V.S. Cheong, B.C. Roberts  fibular proximal growth plate and the tibio-fibular joint, and detecting the largest connected region (MATLAB 2018A, The Math-Works Inc., Natick MA, USA). For cases where the automatic removal of the fibula is incomplete, a mask was created to aid the script in the removal of the fibula. After alignment, the images were cropped to 80% of the tibial length starting from the slice adjacent to the most distal section of the proximal growth plate, to exclude changes in the growth plate in the analysis [26] . Density calibration constants provided by the micro-CT manufacturer, with weekly quality checks using a five-rod densitometric phantom (Scanco Medical, Brüttisellen, Switzerland), were used to convert the linear attenuation coefficients at each voxel to tissue mineral density (TMD) values. Thereafter, the micro-CT images were segmented by applying a single level threshold to binarise the images. This was determined as the midpoint between the background and bone peaks in the grey value frequency plot [25] .
The threshold values used in this study were 498.3 ± 38.6 mg/cc. The subject-specific threshold value computed at week 14 was applied throughout the study, to ensure consistency when comparing between predicted bone adaptation (pseudo micro-CT images) and the experimental dataset ( Section 2.5 ).

Micro-FE analysis
Micro-FE models for the micro-CT images acquired at weeks 14, 18 and 20 were created by directly converting all bone voxels from the cropped images into linear hexahedral elements (voxel-based brick elements). This mesh generation technique is computationally efficient, and directly links the microstructural properties from the micro-FE models to the observed bone remodelling changes. All elements were assigned homogenous isotropic material properties ( E = 14.8 GPa, ν = 0.3), which has been validated to reproduce the stiffness and displacement that were experimentally measured using in situ mechanical testing and DVC [28] . The boundary conditions set in the micro-FE models were used to approximate physiological loading ( i.e. components of loads in the axial and posterior-anterior directions, which induced bending and compression due to the curvature of the tibia [29] ), to investigate the changes in the micromechanical properties due to OVX and after treatment. The model was fully constrained at the proximal end, while the distal nodes were connected with kinematic coupling to the area centroid of the distal surface to avoid unrealistic bending moments during the applications of the axial and posterior-anterior loads ( Fig. 1 B) [11] . A 1 N static load was applied to the control node separately in the axial, posterior-anterior and lateral-medial directions for each mouse model for computational efficiency, which was solved on the University of Sheffield high performance computing clusters (ShARC) using Abaqus 2017 (Dassault Systèmes Simulia, USA). All required input files for the micro-FE models were generated with Matlab and each model contained approximately 9 million elements. The results were postprocessed through scaling and superimposition of the effects to peak physiological loading [30] , according to the body mass (BW) of the mouse at each time point, as described in Cheong et al. [11] (0.01355 * BW N/g along the superior-inferior direction and 0.00289 * BW N/g along the posterior-anterior direction; a coefficient of 0 was applied to the medial-lateral direction as it has minimal effect on the strain energy density). The average strains along the tibial length were plotted to determine the frequency distribution and changes in trend.

Computational algorithm
Simulations of bone adaptation were conducted by applying a linear mechanoregulation algorithm to the micro-FE models [11] , to determine the extent the changes in structural properties of the  mouse tibia is linearly strain driven. The bone adaptation algorithm is based on the mechanostat theory that bone adapts to changes in mechanical stimuli by altering its TMD [31] . Strain energy density (SED) was used as the driving mechanical stimulus as it has been shown to give realistic result in predicting bone changes [ 2-4 , 32 ]. Changes in TMD were computed by comparing the local stimuli at each node with the remodelling law ( Fig. 2 A), which were applied to background and bone surface voxels respectively. Voxels with an

ARTICLE IN PRESS
JID: ACTBIO [m5G; 9:23 ] initial TMD value around the threshold transition zone (2-3 layers from the surface) were permitted to undergo both apposition and resorption, to account for potential errors in segmentation due to partial volume effect. A total of three parameters in the mechanoregulation algorithm (SED threshold k , rate for bone apposition B app , and rate for bone resorption B res ) can be optimised for each simulation to determine the effects of mechanical loading on the bone regulatory process. The application of mechanical loading can increase the structural stiffness of bone [7] , hence the parameters of the algorithm were optimised to best match the changes in the volumetric second moment [11] . In this paper, the models generated from the micro-CT images acquired at weeks 14, 18 and 20 were used to predict the voxel mesh at weeks 16, 20 and 22 ( in silico dataset), respectively, which were compared with the voxels mesh segmented from the micro-CT images acquired at weeks 16, 20 and 22 ( in vivo dataset). A pure mechanoregulation algorithm was used to predict both the effects of OVX (week 14-16) and mechanical loading (weeks 18-20, 20-22) as it was found to be able to predict densitometric changes due to OVX with no significant difference between the experimental and predicted results in the murine caudal vertebra [ 3 , 4 ].

Evaluation metrics
Surface bone adaptation (remodelling) that occurred between consecutive sets of micro-CT scans were computed by comparing the changes within the endosteal and periosteal outlines of the segmented images, and totalling the number of voxels (frequency) between the two sets of scans that had undergone increases (apposition) or decreases (resorption) in TMD across 10 sections of the bone. To account for the change in tibial length with age due to growth, the images of two consecutive scans were superimposed by aligning their volume centroids, and cropping the images to obtain the same dimensions.
Comparisons were made between weeks 14-16 and 16-18 to evaluate OVX effects, and between weeks 18-20 and 20-22 to evaluate load-induced changes. Densitometric parameters bone volume (BV), bone volume fraction (BV/TV), bone mineral content (BMC) and volumetric bone mineral density (BMD) were computed for both the in vivo and in silico dataset, as detailed in Lu et al. [27] . These analyses were conducted for the whole bone and 40 compartments (10 longitudinal sections, 4 partitions in each section in the anterior, posterior, lateral and medial regions) ( Fig. 2 B).
To compare the accuracy of the locations of the predicted bone adaptation sites, which can be interpreted as the extent the bone changes are linearly strain driven, three metrics were computed [11] : 1) the overlap ratio as the volume of intersection between the two datasets after binarisation, normalised by the total volume of the two datasets; 2) the spatial match as the correct bone changes in apposition and resorption in each voxel (a measure of type I error); 3) the prediction accuracy as the number of voxels predicted correctly by the model, normalised by the experimentally measured bone changes (a measure of type II error) ( Fig. 2 C-E). The locations of surface bone changes were also visualised using 2D and 3D images to compare the results qualitatively.

Statistical analysis
The Wilcoxon signed-rank test was used to test for any significant difference between two time points within each dataset, while the Mann-Whitney U test was conducted to compare the effect of treatment at each time point between the treatment group and the control groups (Origin 2018, OriginLab Corp., Northampton, MA). Heat maps representing the p values are presented across the 10 sections of the analysis tibia and statistical significance was set at p < 0.05 (two-tailed).

Densitometric changes due to ovx and mechanical loading
To separate the effects due to OVX, and the impact of mechanical loading on top of OVX, the data from the longitudinal analyses have been separated into two parts and summarised in Figs. 3 and 4 , respectively. BMC increased in WT mice from weeks 14-18, but had modest changes in OVX and ML mice (both ovariectomized), and caused the significant difference in BMC between WT and ML groups at week 14 to become insignificant at weeks 16 and 18 ( Fig. 3 A). A similar trend was observed in BMD at weeks 14 and 16. At week 18, there was however a significant difference between WT and OVX mice (effect of OVX), and unexpectedly between OVX and ML mice (highlighting the importance of normalising the data before the treatment starts).
The external mechanical load applied at week 19 triggered significant increases of 16.1 ± 4.0% and 11. 5 ± 4.0% in BMC and BMD between weeks 18-20, respectively (OVX vs ML, Fig. 4 A). Mechanical loading at week 21 also increased BMC and BMD from weeks 20-22, but only the 6.0 ± 3.0% increase in BMC was significantly different. The first week of treatment countered the effects of OVX as there was a significant difference in BMD between the ML and OVX groups at week 20, which was tending towards significance for BMC ( p = 0.052). The second applied load at week 21 also countered the effects of OVX despite the pause in loading at week 20, as BMC was significantly different, while BMD was tending towards significance between the two groups ( p = 0.052).

Bone adaptation due to ovx and mechanical loading
The most active region of bone adaptation occurred in the most proximal region (section 10) of both endosteal and periosteal surfaces across all time points and groups ( Figs. 3 B and 4 B). OVX increased endosteal resorption more than it decreased endosteal apposition in both weeks 14-16 and 16-18. Periosteal apposition was higher in the OVX and ML groups than the WT group at both time points. Periosteal resorption was also higher in the OVX-operated mice than the WT, with more longitudinal sections showing significant differences at weeks 16-18 than weeks 14-16.
In the ML group, periosteal resorption was approximately 40% higher in the diaphysis (sections 4 to 7) at weeks 16-18 than at weeks 14-16. Endosteal apposition was also higher during the later period of OVX with the highest change of approximately 60% in the distal regions (sections 1 and 2) and the diaphysis (sections 6 and 7).
After the application of passive mechanical loading at week 19, there was a 90-150% increase in endosteal apposition in the bone diaphysis (sections 4 and 5), and a 50% reduction in endosteal resorption in the proximal tibia (sections 1 to 3) at weeks 18-20 than that at weeks 14-16. There was a downward shift in periosteal resorption, affecting almost the entire tibial length, with the highest reduction of approximately 60% in the diaphysis (sections 5-7). Periosteal apposition at weeks 18-20 was more localised in the diaphysis which increased by 50-80% (sections 5-8). Between weeks 16-18 and 18-20, and compared with the OVX control, the effect of mechanical loading was higher on the periosteal section than on the endosteal section.
At weeks 20-22, after the second application of mechanical loading at week 21, endosteal resorption was higher than that at weeks 14-16 between the diaphysis and the proximal metaphysis (sections 7-9), but lower in the diaphysis (sections 4 and 5). Endosteal apposition at weeks 20-22 was significantly lower than that at weeks 18-20 across most of the tibial length.   8-10), which were up to 30% lower than those at weeks 14-16. Periosteal resorption at weeks 20-22 remained lower than that at week 14-16, but was up to 70% higher than that at weeks 18-20 in section 8. Similarly, periosteal apposition remained lower than that at weeks 14-16, but was slightly lower than that at weeks 18-20 by approximately 10%. There were fewer longitudinal sections with significantly higher apposition or lower resorption between the ML and OVX in the first and second sets of mechanical loading, than between pre-loading and after the first set of loading.

Strain distribution predicted by the FE models
The analyses of the SED values in the models at the different time points ( Figs. 3 C and 4 C) showed the effects of accelerated bone resorption and bone formation due to OVX and subsequent mechanical loading of the mouse tibia, respectively, leading to changes in structural properties. Average SED generally increased as a result of OVX from weeks 14 to 18. The difference in SED distribution between OVX and WT mice were significantly different in most longitudinal sections of the tibia at week 14, but became largely similar at weeks 16 and 18 ( Fig. 3 C). The applied loads at weeks 19 and 21 reduced the average SED across the tibia, which was significantly different to the OVX groups across all 10 longitudinal sections ( p < 0.01). There were more longitudinal sections associated with a significant reduction in SED after the first set of mechanical loading (ML1), than after the second set of loading (ML2). The reduction in SED was also slightly lower after the second set of mechanical loading, with the exception of the proximal metaphysis (90-98% of tibial length; section 10), which also had the largest reduction in SED (40%). Persistent reduction in SED after each set of loading was observed only in the distal region, and sections 6 and 9.
The strain distribution in the tibia of a representative mouse, for normalised physiological loads (by its body weight BW), showed regional differences in the values of SED, maximum principal strain (tensile) and minimum principal strain (compression) across the bone ( Fig. 5 ). The highest magnitudes of maximum principal strain were found near the interosseous crest (IC) in the distal regions, and between the proximal tibial crest (PTC) and tibial ridge (TR) in the proximal regions. The highest magnitudes of minimum principal strain were located around the IC at the proximal region, and at the medial region distally. The SED results, accounting for the different components of the strain, as expected showed the highest magnitudes in all the sections described above. From weeks 14 to 16, high increases in SED were located in the TR of the proximal metaphysis (section 9-10) and the medial distal region (5-15%) ( Figs. 5 and 6 ). After the first application of mechanical loading, high SED reductions were observed on the medial aspect of the distal tibia and at the anterior aspect and IC of the proximal tibia, which increased slightly at week 22 ( Fig. 5 ).
SED frequency plots of 5 longitudinal sections showed that the distribution is region-specific ( Fig. 6 ). The proximal regions (sections 7 and 9) also displayed a wider range of SED values compared to the diaphysis (sections 3 and 5). Bone adaptation induced by the applied load caused the distribution to become narrower with a shift towards lower SED values. This 'pivoting' in the graphs differed for each region. For example, the intercept between the graphs before loading at week 18

Bone adaptation parameters
The parameters assigned to the model showed that the first application of mechanical loading at week 19 significantly in-creased the apposition rate ( p < 0.05), and reduced the remodelling threshold ( p < 0.05), compared to the parameters computed between weeks 14-16. A trend in the reduction of the resorption rate was also observed, but this was not significantly different ( Table 1 ). The same trend was also observed after the second application of mechanical loading at week 21, but the change was smaller and not significantly different to the values found between weeks 18-20. All the parameters of the model obtained after 2 applications of mechanical loading (between weeks 20-22) were significantly different to those calculated between week 14-16, when it was under the effect of OVX alone. ( p < 0.05).

Simulation accuracy
The in silico models were able to capture the trends shown by the densitometric measurements. The predicted BV at weeks 16 and 20 were 9.4 ± 0.4 mm 3 and 10.4 ± 0.4 mm 3 , respectively, which were not significantly different from the experimentally measured values of 9.4 ± 0.5 and 10.5 ± 0.5 mm 3 , respectively. Similarly, the predicted BV/TV at weeks 16 and 20 were 0.56 ± 0.01 and 0.61 ± 0.01, not significantly different from the experimentally measured values of 0.56 ± 0.01 and 0.59 ± 0.01. Both BV and BV/TV were over predicted at week 22 ( p ≤ 0.016). The BV at week 22 was over predicted by 3.8% from 10.9 ± 0.5 mm 3 , while the BV/TV was over predicted by 5.2% from 0.60 ± 0.02.
The predicted BMC values at weeks 16 and 22 were 9.4 ± 0.4 mg and 11. 9 ± 0.7 mg, respectively, which were not significantly different from the experimentally measured values of 9.4 ± 0.5 mg and 11. 5 ± 0.4 mg, respectively. The predicted BMD values at weeks 16 and 20 were also not significantly different from the experimentally measured values. The BMC at week 20 was under predicted by 4.4% from 10.9 ± 0.5 mg and the BMD at week 22 was over predicted by 4.3% from 0.64 ± 0.02.
The distributions of the errors for the different partitions were similar amongst the densitometric parameters, thus only the results for BMC and BMD will be described in this section ( Fig. 8 ). The mean absolute errors were below 20% across the 10 sections (including all 4 compartments, Fig. 8 E, J), with the lowest errors at week 16 ( < 10% error), within the cortical bone (sections 1-9). Across all 4 compartments at week 16, the highest errors were produced at the proximal end with trabecular bone (section 9 and 10), which were up to 15.7% for BMC and 13.9% for BMD. At week 20, errors of up to 21.3% were also produced at the proximal end across all 4 compartments for BMD. For BMC, the proximal regions showed the highest errors (26.3%) in the medial section. The overall errors reduced at week 22, with the highest errors of 29.5% in the proximal medial compartment for BMC, and 24.3% in the proximal anterior compartment for BMD.
The lowest BMC errors at weeks 16, 20 and 22 were 0.2% in section 6 of the medial compartment, 1.0% in section 8 of the posterior compartment and 0.5% in section 6 of the posterior compartment, respectively. The lowest BMD errors at weeks 16, 20 and 22 were 0.4% in section 6 of the medial compartment, 0.1% in section 4 of the lateral compartment and 0.1% in section 5 of the anterior compartment. Similar trends were found for BV and BV/TV (Fig.  S1).
Visual analysis of the evaluation metrics at 30%, 50%, 70% and 90% (sections 3, 5, 7 and 9) of the analysed tibia length for a representative mouse ( Fig. 9 ) showed high degree of overlap between the predicted and experimental images for apposition. The spatial match for all three models were higher on the periosteal surface than on the endosteal surface for apposition, but the opposite was true for resorption. The spatial match and predictive accuracy across the mouse tibia displayed different patterns at all three time points on both the endosteal and periosteal surfaces (Fig. S2)

Discussion
The goal of this study was to measure the effects of OVX and subsequent mechanical loading treatments on bone adaptation in the mouse tibia, and to quantify how much the load-induced bone changes are correlated with changes in the tibial geometry. This is the first time the location and extent of bone adaptations are quantified throughout the tibial length after two sets of mechanical loading, in an OVX tibial model. Moreover, micro-FE analyses were Please cite this article as: V.S. Cheong, B.C. Roberts   conducted at each time point to determine the local strain changes that resulted from the progression of OVX and the effect of the applied load, to link the contribution of strain changes to bone adaptation and modifications in the tibial shape. This study also benefits from the modelling of bone adaptation, and the use of the volumetric second moment as the optimisation criteria between consecutive set of images in the bone adaptation algorithm, to determine the extent the geometrical changes were linearly strain driven.
This study uses a simple loading model under physiological loading coupled with a bone remodelling algorithm to determine the effects of the interventions on the local tissue changes. The spatial match (week 18-20: 69 ± 4%, week 20-22: 61 ± 2%) and prediction accuracy (week 18-20: 53 ± 6%, week 20-22: 44 ± 3%) obtained in this study are comparable to the results of Pereira et al.
who obtained a Kendall τ rank coefficient of 0.51, and had applied the load used during the experimental session to the joints [6] . The accuracy of this modelling approach is also comparable to the results of Schulte et al. who obtained an overall spatial match of 55% and 48% in the non-OVX loaded, and the OVX caudal vertebrae, respectively [4] .
The results showed regional differences in bone adaptation on the endosteal and periosteal surfaces during OVX and after each application of mechanical loading in the mouse tibia ( Figs. 3 and  4 ). During the progression of OVX, there was an increase in periosteal resorption in weeks 16-18 compared to weeks 14-16, while the curves for endosteal remodelling and periosteal apposition were similar. In contrast, healthy mouse tibia between weeks 14 to 22 exhibited higher periosteal apposition, and a slight reduction in periosteal resorption in the proximal metaphysis with age [11] .
Densitometric parameters following OVX showed an initial reduction from weeks 14 to 16, before a slight increase between weeks 16 to 18, although the differences were not significant ( Fig. 3 ). Previous longitudinal study using the OVX caudal vertebra model showed an initial increase in BV/TV before reducing with time [ 4 , 15 ]. However, a later study showed an initial decline, before a small increase in BV/TV 9 weeks after OVX [3] , similar to the results obtained in this study.
Following OVX, morphological changes caused by the overall net increase in periosteal resorption than endosteal apposition, and the decreases in BV and BV/TV, translated to higher average strains under physiological loading ( Fig. 3 , Figs. 5 and 6 ). The increase in SED was higher between weeks 14 to 16 than weeks 16 to 18, and as the bone changes predicted by the combination of homogeneous micro-FE models and mechanoregulation algorithms are driven by geometrical changes [ 4 , 7 ], the results indicate that most of the structural changes occurred early on. The increase and decrease in endosteal apposition and resorption, respectively, in the diaphysis and the proximal regions (sections 6-8), could have resulted in the slight increase in densitometric parameters between week 16 and 18. Moreover, as net bone formation on the endosteal surface is less effective than bone formation on the periosteal surface in resisting bending, this would explain the similarity in the average SED in the diaphysis (sections 5-8). The slightly lower densitometric properties and higher strain distribution in the OVX-operated mice than the WT group in this study are consistent with previous studies that compared the tibiae of OVX and SHAM mice [18] .
Mechanical loading at week 19 caused increases in both endosteal and periosteal apposition, and reduction in both endosteal and periosteal resorption, in line with recent findings in the literature [ 8 , 12 ]. The effects of mechanical loading at week 21 on bone remodelling were much reduced, as the periosteal resorption was higher, while periosteal apposition was slightly lower than at weeks 18-20 ( Fig. 4 ). There was little impact on endosteal remodelling, as endosteal apposition returned to the levels observed between weeks 14-16, and endosteal remodelling was higher in some regions than that at weeks 14-16. In contrast for healthy mice, mechanical loading of the tibia at 12 weeks, over a 2-week period, showed a reduction in the rate of bone formation between weeks 13 and 14, compared to weeks 12 and 13 [12] . Interestingly, the eroded surface was found to be lower between weeks 12 and 14 than between weeks 12 and 13. The observed differences in results are probably due to the differences in the time of application of the mechanical loading. In fact, Javaheri et al. [12] applied the loads at weeks 12 and 13 when the tibia is still growing [26] , while in this study the loads were applied at weeks 19 and 21 of age. Visualisation of the locations of bone adaptation showed that bone apposition in the proximal tibia occurred on both endosteal and periosteal surfaces after the first period of mechanical loading between weeks 18 and 20 ( Fig. 9 ). Between weeks 20 and 22, it was found that apposition occurred on opposite surfaces (medial periosteal and lateral endosteal). The amount of bone apposition in the proximal tibia after one period of loading is similar to the results obtained by Javaheri et al. [12] , who applied a similar loading rate in 12-weeks-old mice and found bone apposition throughout the proximal tibia after two weeks of mechanical loading, and larger increases of bone formation in the first week than the second. The result in this study is higher than the amount of apposition reported in the tibia of healthy 26-weeks-old mice, as apposition was shown to be located primarily on the posterior endosteal surface in the proximal tibia, and on the medial and lateral endosteal surfaces at 50% of the tibial length [8] . The differences could be due to the age of the mice used in the experiments and the application of the mechanical loading on ovariectomised or non-operated mice.
Bone resorption was low between weeks 18 and 20, which was similar to the results obtained by Birkhold et al. [8] for 26-weeksold mice. However, increases in bone resorption between weeks 20 and 22 occurred on the medial endosteal surface at the proximal end (70-90% of tibial length), and the medial surfaces and lateral endosteal surface at the distal end (30-50% of tibial length), but in young healthy mice this was found to be located primarily on the posterior aspect [12] . In particular, the pattern obtained at 50% of the tibial length after the second application of mechanical loading, between weeks 20 and 22, was similar to the pattern obtained by Javaheri et al. [12] .
The overall outcome is that loading had a slightly greater effect in reducing bone resorption, than on increasing bone formation, similar to the results reported in peri-implant bone formation in the caudal vertebrae [33] . This could be due to the higher overall strain in the bone after the first set of loading, which outweighed the signals for systemic bone resorption, leading to a reduction in resorption after the first week of loading. As the bone becomes stiffer, the surface strains induced by the week 21 load decreases, leading to lower signals to counteract systemic signals for bone resorption and for targeted bone formation.
The changes in bone adaptation due to mechanical loading ( Fig. 3 ) translated to the significant increases in densitometric parameters ( Fig. 7 ), similar to previous studies that demonstrated the anabolic effect of mechanical loading in OVX bones [ 17 , 18 ]. In the whole OVX caudal vertebra, the second mechanical loading led to a higher rate of increase in BV/TV than the first [3] . The discrepancy could be due to the greater amount of trabecular bone present in the caudal vertebra than in the mouse tibia. In contrast, a smaller increase in BV/TV after the first loading in healthy mice has also been shown [4] . These morphological responses to mechanical loading, in particular cortical thickening [ 8 , 12 , 16 ], were also reflected through the reductions in strain distributions from weeks 18 to 22 ( Fig. 4 -6 ). This was found to be more dominant in the proximal tibia during the first period of loading and in the diaphysis during the second period of loading [12] . However, the SED distribution from this study ( Fig. 6 ) showed that the highest reductions in strains were in the distal tibia (0-20%) and the prox-imal tibia (80-100%) after the first load-induced gain, whereas reduction in strain occurred mainly in the proximal tibia (85-98%) after the second mechanical loading. Previous results have shown that increases in BMC and cortical thickness induced by prior mechanical loading were not maintained 4 weeks after OVX [18] . The results from this study show that skeletal loading 4 weeks after OVX were still able to achieve gains in densitometric and structural properties to compensate and exceed the loss sustained during OVX. The optimised bone adaptation parameters showed a reduction in remodelling threshold, bone resorption rate, and an increase in bone apposition rate after each application of mechanical loading ( Table 1 ), to predict the overall increase in bone apposition observed experimentally between weeks 18 and 22 ( Fig. 3 ). The pa-

ARTICLE IN PRESS
JID: ACTBIO [m5G; 9:23 ] rameters obtained in this study are different from those obtained in the caudal vertebra [ 3 , 34 ] as they had used the same remodelling threshold and remodelling rate for both periods of mechanical loading, by optimising the predictions to the average response of BV/TV. The in silico models have been able to capture the same trends in the densitometric measures ( Fig. 7 ) and the predictive accuracy of the model was over 40% in OVX and over 60% after mechanical loading ( Fig. 8 ). The errors in predicted densitometric measures are similar to the errors reported by Levchuk et al. [3] and Schulte et al. [4] . Hence the results suggest that the reduction in measured bone resorption at week 20 was captured primarily by the reduction in the remodelling threshold. Moreover, as the increase in the volumetric second moment (cortical thickness) caused overall strain measures to decrease between weeks 20 and 22, while bone apposition remained high, this then caused the remodelling threshold and the resorption rate to decrease from weeks 20 to 22 in compensation for the reduction in strain. The high spatial match and predictive accuracy in apposition were not matched in resorption. Moreover, the models over predicted densitometric measures after mechanical loading, despite optimising the parameters every 2 weeks to capture any potential change in trend. The over prediction was higher at week 22 than at week 20, which could be due to the stronger combined effect of age and OVX leading to higher periosteal resorption, and which the model had not accounted for. These results show that bone apposition is strain driven at the organ-level, but not bone resorption. The regional differences in the SED distribution along the tibial length suggest that the threshold for remodelling may be section dependant, and hence resorption could be more locally driven, as locations of resorption have been found to be more sporadically located [19] . Moreover, a recent study using deep neural network to assess the rejuvenation effect of repeated mechanical loading in a two-week period showed that the greyscale values in the region immediately below the growth plate were more important for estimating the age of the bone at day 0 and after four repeated sessions of mechanical loading, but shifted to the distal portion of the proximal tibia (corresponding approximately to section 8 in this paper) after 7 and 10 sessions of loading [35] . The higher accuracy in predicting load-induced apposition than during OVX showed that most of the bone response to mechanical loading is indeed linear, but up to 30% of the changes in the cortical bone were independent of strain. Non-strain response was higher at 50% in the proximal tibia, where trabecular bone was present. As the adaptive capacity of bone in response to load decreases with age [19] , the overall predictive accuracy of the model in predicting load-induced apposition was higher at week 22 than at 20. All these reasons might explain the frequent over-prediction of bone ingrowth in implants using the mechanoregulation models [ 2 , 32 ]. Future work could include a long detraining period to investigate the frequency of exercise required, and predicting the bone adaptation over a longer period, to make the findings of animal studies more translatable to humans.
A small sample size ( N = 6) was used, which is consistent with other longitudinal imaging studies that measured bone changes in individual mice over time [33] . As changes are normalised to baseline values to account for any difference prior to the start of the analysis, longitudinal studies benefit from lower variance compared with cross-sectional studies, and have been able to demonstrate significant difference even with a small sample size [ 26 , 36 , 37 ]. Indeed, the BMC and BMD values at baseline were significantly lower for the WT group than the ML group used in this study, and there was also significant difference in the BMD between the OVX and ML group at week 18. The differences could in part be due to different generations of mice, hence most of the values reported are changes relative to the baseline.
There are several limitations in this study. Firstly, there was some variability in the adaptive response to mechanical loading, which may be due to the positioning of the tibia in the loading jig during in vivo loading [10] , as the fixture uses a semi-spherical cup in line with other studies [20] . Moreover, the strains applied at week 21 were not strain matched to those at week 19, which may have resulted in a lower adaptive response to mechanical loading. Nevertheless, differences in the surface remodelling response and the strain distribution between the 2 sets of loading in Fig. 4 show that adaptive changes were more localised after the second set of loading. Hence, applying a higher load that is strain-matched at week 21 may not result in a similar adaptive response as between weeks 20-22. The experimental variations have also not been accounted for in the micro-FE model, making it difficult to distinguish if the variations in the remodelling parameters were due to experimental variability or differences in response to OVX and/or mechanical loading. Validation of the induced strain in vivo has not been conducted in this study [20] , as the application of strain gauges may cause a local stiffening of the specimen, as shown for the mouse forearm [38] , and their application may affect the success of the longitudinal study due to artefacts in the in vivo micro-CT images, increased risk of mortality due to local sepsis etc. Nevertheless, longitudinal monitoring of the murine tibiae throughout this study allows for comparison to be made between the micro-FE models to evaluate the effects of OVX and/or treatment, as the models were simulated with the same boundary conditions and time-point adjusted loads. Voxel-based brick elements were used for the micro-FE models, which does not capture the surface strain as well as tetrahedral elements. However, this was mitigated with the implementation of a bone remodelling unit lattice which included the strains of 2-3 layers of voxels from the surface, coupled with an optimisation algorithm to find the parameters that provided the closest match to the experimental results [11] . A nodebased approach was utilised to prevent checkboard problems of discontinuities, but it adds an additional step to extrapolate the results calculated at the gauss points to the nodes during postprocessing [39] . The applied mechanoregulation algorithm assumes a linear response of bone adaptation to local strain changes, but there could be other non-linear changes that this study has not explored. Furthermore, this algorithm cannot distinguish between modelling and remodelling, and the ability to separate the two may improve the predictive accuracy of the model. The loads and boundary conditions were applied to the cropped bone rather than through the joints, as the material properties of the tibio-fibular growth plate and joint are not known [11] , and the effect of the removal of muscles and soft tissues of the lower leg on load transmission is small [40] . Hence, the FE models were used in this study to compare the changes in the trend due to disease and mechanical loading, and not the absolute strain distribution. Future work should include sensitivity analysis on the effect of the presence of the growth plate and the fibula on the strain [41] . Homogenous micro-FE models were used throughout the study as the focus of this study was on bone morphological changes. Despite not including potential changes in local bone mineralisation, micro-FE models with homogenous material properties have been shown to predict the experimentally measured stiffness and local displacement values accurately [28] . Moreover, the differences between the average strain distribution obtained from heterogeneous and homogenous micro-FE models have previously been reported to be minimal across the tibial length in adult mice [7] . However, the use of heterogeneous material properties may capture local strain in the newly formed bone due to mechanical loading better than homogenous models. This approach will be tested in the future to evaluate potential improvement in the prediction of local SED and in bone remodelling. The use of deep learning approaches has demonstrated high accuracy in predicting skeletal age from micro-Please cite this article as: V.S. Cheong, B.C. Roberts and V. Kadirkamanathan et al., Bone remodelling in the mouse tibia is spatiotemporally modulated by oestrogen deficiency and external mechanical loading: A combined in vivo/in silico study, Acta Biomaterialia, https://doi.org/10.1016/j.actbio.2020.09.011 CT images in both loaded and unloaded control, and hence could potentially be incorporated in this model to separate the effects due to modelling and remodelling better [35] . The amount of resorption predicted in this study is low, and while it is similar to previous studies that reported on the underestimation of bone resorption [ 3 , 4 ], it shows that this algorithm is unable to predict bone adaptation driven by biological changes. Future work should include the development of a multi-scale framework to account for biochemical signals explicitly using biological networks or agentbased models [42] . The comparison of the longitudinal images depends highly on the accurate registration of each set of micro-CT scans, and a validated registration approach that produces errors of less than 3.5% was used in this study [ 26 , 27 ]. However, some of the results ( e.g. between weeks 20-22 in Fig. 9 ) suggest that registration errors may be more concentrated in certain regions than others. Finally, the higher bone turnover in the trabecular bone led to poor prediction of densitometric measurements in the proximal tibia. This suggests the need to have two or more separate models, especially prior to week 18, as the remodelling threshold seems to be region dependant.

Conclusion
Using in vivo longitudinal micro-CT imaging to monitor the effects of mechanical loading treatment in OVX tibia, revealed for the first time that the load-induced response was higher during the first set of mechanical loading, and also showed an increase in periosteal resorption after the second set of mechanical loading. The load-adaptive response led to a decrease in average strain across the mouse tibia after the first loading, but the strain distributions in the tibia were similar between the first and second loading. This is the first bone adaptation simulation of mechanical loading in the mouse tibia with OVX, and the results showed a linear bone apposition response to load-induced changes, which was higher during mechanical loading than during OVX. The results also showed that resorption was independent of strain at the organ level, but a region-dependant distribution of strain suggests that resorption may be locally modulated or biologically driven, or both. A higher loading rate or peak load may thus be required to modulate the increase in resorption with age.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.