Mechanostat parameters estimated from time-lapsed in vivo micro-computed tomography data of mechanically driven bone adaptation are logarithmically dependent on loading frequency

Mechanical loading is a key factor governing bone adaptation. Both preclinical and clinical studies have demonstrated its effects on bone tissue, which were also notably predicted in the mechanostat theory. Indeed, existing methods to quantify bone mechanoregulation have successfully associated the frequency of (re)modeling events with local mechanical signals, combining time-lapsed in vivo micro-computed tomography (micro-CT) imaging and micro-finite element (micro-FE) analysis. However, a correlation between the local surface velocity of (re)modeling events and mechanical signals has not been shown. As many degenerative bone diseases have also been linked to impaired bone (re)modeling, this relationship could provide an advantage in detecting the effects of such conditions and advance our understanding of the underlying mechanisms. Therefore, in this study, we introduce a novel method to estimate (re)modeling velocity curves from time-lapsed in vivo mouse caudal vertebrae data under static and cyclic mechanical loading. These curves can be fitted with piecewise linear functions as proposed in the mechanostat theory. Accordingly, new (re)modeling parameters can be derived from such data, including formation saturation levels, resorption velocity moduli, and (re)modeling thresholds. Our results revealed that the norm of the gradient of strain energy density yielded the highest accuracy in quantifying mechanoregulation data using micro-finite element analysis with homogeneous material properties, while effective strain was the best predictor for micro-finite element analysis with heterogeneous material properties. Furthermore, (re)modeling velocity curves could be accurately described with piecewise linear and hyperbola functions (root mean square error below 0.2 µm/day for weekly analysis), and several (re)modeling parameters determined from these curves followed a logarithmic relationship with loading frequency. Crucially, (re)modeling velocity curves and derived parameters could detect differences in mechanically driven bone adaptation, which complemented previous results showing a logarithmic relationship between loading frequency and net change in bone volume fraction over 4 weeks. Together, we expect this data to support the calibration of in silico models of bone adaptation and the characterization of the effects of mechanical loading and pharmaceutical treatment interventions in vivo.

per running execution). Since the pipeline is faster, a higher number of samples can be processed in a very short period. Moreover, as each sample is analyzed in a single core, our parallel implementation further reduces the processing time, being only limited by the number of cores available. For instance, the images from a 4-week study comprising 40 samples per time-point would require four iterations on a 48-core CPU to perform this compartmental analysis, in less than 10 minutes, or the time equivalent for the IPL version to analyze a single sample. Even if the IPL implementation would process ten images in parallel, the proposed Python version would still represent at least a 16x speed up in processing time. More importantly, the proposed implementation can run on the desired HPC platform, ideally where subsequent micro-FE and post-processing steps will also be performed, reducing auxiliary data transfer operations to move data between systems, which would further increase the total overhead of the analysis and disk space requirements. Additionally, as the outputs generated for each sample are exported to a single file, the data management required to oversee all processed samples is considerably simplified compared to the IPL version.

Micro-finite element analysis pipeline
As a second step, the micro-FE pipeline implemented in Python was also validated against the established method in IPL. Prior to the micro-FE, simulated intervertebral discs must be added at the top and bottom of the vertebra (Webster et al. 2008). The output of the proposed and established implementations achieved an average DSC of 0.997 (median: 0.998, interquartile range: 0.996 -1.000). Furthermore, regarding the execution time for the default settings for all samples (stiffness: 14.8 GPa and radius: 13% of the height of the vertebra), the IPL implementation needed an average of 15.81 seconds (median: 16 sec., interquartile range: 15 sec. -17 sec.) of CPU time per sample to add the intervertebral disc to the image, running on HP OpenVMS (version V8.4) with HP Integrity BL890c i2 (1.60 GHz). The proposed Python method required an average of 28.18 seconds (median: 27.60 sec., interquartile range: 26.64 sec. -28.85 sec.) of CPU time per sample running on Fedora 35 with Intel Xeon Platinum 8168 CPU (2.7-3.7 GHz). Despite the slight increase in execution time, there are key advantages of the Python implementation, mainly that while the IPL version is designed to work exclusively on binary images, therefore only supporting its use for micro-FE analysis using homogeneous material properties, the proposed Python version expands this approach to grayscale images.

Calibration of simulated intervertebral disc parameters
A preliminary analysis assessed the effect of the Young's modulus value and radius of the simulated intervertebral discs added to the vertebra on the ability to quantify mechanoregulation information from time-lapsed micro-CT data and identify the most suitable parameters to be used for the micro-FE analysis with homogeneous and heterogeneous material properties. The default radius on the established IPL pipeline is 13% of the height of the vertebra (Webster et al. 2008), and previous work has used stiffness values of 2.5 MPa (Webster et al. 2008) and 14.8 GPa (Schulte et al. 2013). Therefore, a disc radius between 7% and 17% of the height of the vertebra and ten stiffness values regularly spaced on a logarithmic scale between 2.5 MPa and 20 GPa (and including 14.8 GPa) were considered. Given the considerable amount of parameter combinations, only the samples of the static loading group from the first time-point were considered, as this is the in vivo group that best mimics the loading condition of the micro-FE analysis.
A micro-FE analysis was performed for all selected samples, using homogeneous and heterogeneous properties, and the mechanoregulation information was quantified using the correct classification rate (CCR) (Tourolle né Betts et al. 2020) as an indicator of the performance of each.
The best-performing settings for each configuration were used for the mechanoregulation analysis performed in the remaining work.
For the micro-FE analysis with homogeneous material properties, the ability to retrieve mechanoregulation information from time-lapsed data showed slight variation (below 1%) across the domain space. Nonetheless, a stiffness value of 14.8 GPa (the same value assigned to bone) or higher showed slightly higher CCR values. Concerning the disc radius, a value between 11% and 13% of the vertebra height also matched the region of the highest CCR values. Additionally, this setting enabled the solver to converge for all samples in a single run, while lower disc stiffness values required adjusting the solver settings for some samples. Conversely, there was a higher variability of CCR values across the domain for heterogeneous material properties, with a stiffness of 18.4 MPa achieving the highest CCR value, for a radius of 13% of the vertebra height. In both instances, lower radii showed a decreasing trend in the amount of mechanoregulation information recovered. A critical remark is that the relationship applied to convert bone mineral density to Young's modulus yielded structures with an average stiffness ranging from 16.9 ± 4.3 GPa (week 0) to 17.5 ± 4.5 GPa (week 4), which is higher than the corresponding Young's modulus value used for homogeneous material properties.

Supplementary Figures
Supplementary Figure 1. Comparison of bone morphometric parameters computed for the regional compartments generated using the established IPL and the proposed Python implementations. A) Dice similarity coefficients computed between the compartments identified with each method for all samples and all time-points. Data is shown per compartment and based on the relative percentage difference of voxels between compartments, assuming the IPL method as the reference. B) Selected static and dynamic bone morphometric parameters computed for all time-points, using compartments identified with the IPL (dashed line) and proposed Python (full line) methods, showing no statistical differences between both methods (p>0.05), as determined by linear mixed effects model; furthermore, significant differences between groups determined by Tukey post-hoc multiple comparisons test of the values from the last timepoint (p<0.05, a: 2 Hz-5 Hz, b: 2 Hz-Sham, c: 2 Hz-Static, d: 5 Hz-Sham, e: 5 Hz-Static, f: 10 Hz-Sham, g: 10 Hz-Static) are still conserved for the proposed Python implementation. C) Correlation between static morphometry parameters (bone volume fraction -BV/TV, specific bone surface -BS/BV, trabecular thickness -Tb.Th, trabecular spacing -Tb.Sp) and dynamic morphometry parameters (bone formation rate -BFR, mineral apposition rate -MAR, bone resorption rate -BRR, mineral resorption rate -MRR) computed with each method for all samples and all time-points. The Spearman correlation coefficients computed show a near-perfect agreement in morphometric indices.

Supplementary Figure 2.
Sensitivity analysis on the effect of the Young's modulus value used for bone and marrow voxels on the quantification of mechanoregulation information, using micro-FE analysis with homogeneous material properties. Only the static loading group was considered for this analysis. CCR values obtained by SED, effective strain, and ∇SED as local mechanical signal descriptors for the mechanoregulation analysis for weekly intervals. In A), bone tissue was assigned a Young's modulus varying between 2.5 MPa and 20.0 GPa, and marrow voxels (non-bone voxels in the trabecular region) were set to 2 MPa. Specifically, ten equally spaced values on a logarithmic scale between the two extremes were considered for the Young's modulus of bone ( Supplementary Figure 3. A) Comparison between the estimated volume with the proposed mechanostat (re)modeling velocity estimation algorithm and actual volume measured by superimposing consecutive time-points for all (re)modeling clusters identified in one sample, with a 1week interval. Three cases are shown: without any correction to the estimated (re)modeling distance (left), with a correction factor based on the total (re)modeling volume (middle), and with a correction factor specific to each (re)modeling cluster identified (right). The distance transform considered the taxicab metric. For each pair, the Spearman Correlation Coefficient is presented. B) Comparison of the total estimated volume with the proposed mechanostat (re)modeling velocity estimation algorithm and actual volume measured between time-points, for all samples and weekly intervals, including or not the correction factor depicted in A) for a single sample, with the corresponding Spearman correlation coefficient for each event. Without any correction, the (re)modeling distance collected from the (re)modeling surfaces yielded a volume change that did not match the volume of their corresponding clusters nor the total volume of (re)modeling clusters between time-points. Using a global correction, the total (re)modeled volume was accurately captured, but the estimated volume of individual (re)modeling clusters still did not match (A and B, middle), leading to a skewed pairing of remodeling distance and mechanical signal values. Finally, using a cluster-specific correction, the estimation obtained from surface voxels is scaled to reflect the corresponding cluster volume accurately for all samples. The volume of each cluster is computed, as well as the estimation from the surrounding surface voxels. Since the resolution of the micro-CT images does not allow determining with certainty from which voxel a given cluster was formed/resorbed, we assume that all neighboring surface voxels of the cluster contribute equally to that process, which is mathematically translated to a linear scaling of the (re)modeling distance values collected from the surface voxels associated with a given cluster to match the change required to produce the volume of the cluster. C) two-dimensional example of the outcome of the surface estimation performed with the proposed algorithm to compute (re)modeling distances for each surface voxel. Using a distance transform applied on the follow-up and inverse of the follow-up images, the displacement of the surface between time-points can be estimated for formation and resorption clusters to approximate the associated volume change of the (re)modeling cluster. Using the cluster-specific volume correction, we ensure that the predicted volume from the surface voxels matches the computed volume of the cluster. Both images were isolated from a threedimensional micro-CT image and show how formation and resorption clusters look and how the corresponding surface is estimated with our approach.
Supplementary Figure 4. A) Matrices with the p-values for multiple comparisons between week intervals (0-1, 1-2, 2-3, 3-4) of the parameter distributions obtained from the piecewise linear function for each group, estimated using the balanced bootstrapping approach described in the methods (see section "Mechanostat (re)modeling velocity curve and parameter derivation"). Parameter legend (see Materials and Methods and Table 1 for an extended description): Resorption saturation level (RSL), Resorption velocity modulus (RVM), Resorption threshold (RT), Formation threshold (FT), Formation velocity modulus (FVM), Formation saturation level (FSL). Statistical significance was determined with Conover's test corrected for multiple comparisons with a step-down method using Bonferroni-Holm adjustments.
Supplementary Figure 5. A) Matrices with the p-values for multiple comparisons between groups of the parameter distributions obtained from the piecewise linear function for each week interval (0-1, 1-2, 2-3, 3-4), estimated using the balanced bootstrapping approach described in the methods (see section "Mechanostat (re)modeling velocity curve and parameter derivation"). Parameter legend (see Materials and Methods and Table 1 for an extended description): Resorption saturation level (RSL), Resorption velocity modulus (RVM), Resorption threshold (RT), Formation threshold (FT), Formation velocity modulus (FVM), Formation saturation level (FSL). Statistical significance was determined with Conover's test corrected for multiple comparisons with a step-down method using Bonferroni-Holm adjustments.
Supplementary Figure 6. A) Matrices with the p-values for multiple comparisons between week intervals (0-1, 1-2, 2-3, 3-4) of the parameter distributions obtained from the hyperbola function for each group, estimated using the balanced bootstrapping approach described in the methods (see section "Mechanostat (re)modeling velocity curve and parameter derivation"). Parameter legend (see Materials and Methods and Table 1 for an extended description): Resorption saturation level (RSL), (Re)modeling threshold (RmT), (Re)modeling velocity modulus (RmVM), Formation saturation level (FSL). Statistical significance was determined with Conover's test corrected for multiple comparisons with a step-down method using Bonferroni-Holm adjustments.

Supplementary Figure 7.
A) Matrices with the p-values for multiple comparisons between groups of the parameter distributions obtained from the hyperbola function for each week interval (0-1, 1-2, 2-3, 3-4), estimated using the balanced bootstrapping approach described in the methods (see section "Mechanostat (re)modeling velocity curve and parameter derivation"). Parameter legend (see Materials and Methods and Table 1 for an extended description): Resorption saturation level (RSL), (Re)modeling threshold (RmT), (Re)modeling velocity modulus (RmVM), Formation saturation level (FSL). Statistical significance was determined with Conover's test corrected for multiple comparisons with a step-down method using Bonferroni-Holm adjustments.

Supplementary Tables
Supplementary Table 1: Parameters of the mathematical functions fitted to the estimated mechanostat group average (re)modeling velocity curves for all weekly intervals. Data presented as "parameter (IQR)", where the IQR was calculated using the balanced bootstrapping approach described in the methods (see section "Mechanostat (re)modeling velocity curve and parameter derivation") to provide an estimate for the variability of each parameter. Root mean squared error (RMSE) was used to characterize the quality of the fit. The row "Effective strain range" indicates the range of mechanical signal values from which the fit of the mathematical functions was derived. Parameter legend (see Materials and Methods and Table 1