Considerable patient variability in bone structure, strength, and day-to-day external mechanical load poses a severe problem in the clinical assessment and treatment of metabolic bone diseases such as osteoporosis. Diagnosis and bone strength assessment rely heavily on radiographic measures of bone mineral density (BMD). However, sources of error in BMD measurements, i.e., intra- and interpatient variability, make it challenging to attribute measured BMD changes to the actual biological change (Nguyen et al., 1997). Accordingly, the sensitivity and specificity of predicting individual patient’s risk for fracture are low (Trémollieres et al., 2010; Cervinka et al., 2017), especially at the hip where falls play a major role. As a consequence, patients may receive treatment, although only a minority would have suffered from a bone fracture. Although these medications are well-tolerated and safe during large-scale clinical trials, anti-resorptive therapies can result in rare and severe adverse events, including osteonecrosis, hypocalcaemia, and thromboembolism (Chen and Sambrook, 2012). Moreover, current diagnostic approaches fail to identify the specific weak spots in the bone. Therefore, they do not estimate where and how fractures will occur and how a local intervention could prevent them (Schultz and Wolf, 2019).

High-resolution peripheral quantitative computed tomography (HR-pQCT), an emerging diagnostic modality of the peripheral skeleton, allows assessing three-dimensional (3D) bone structure and strength at the trabecular level (MacNeil and Boyd, 2007; Melton et al., 2007; Boutroy et al., 2008; Kazakia et al., 2008; Burghardt et al., 2010; Seeman et al., 2010; MacDonald et al., 2011). More recently, complementary methods have been proposed to computationally monitor 3D bone microstructure changes over time (time-lapse) and calculate local mechanical loading using micro-finite element (micro-FE) analysis. This has been demonstrated in mice (Schulte et al., 2013; né Betts et al., 2020; Malhotra et al., 2021) and patients (Christen et al., 2014; Mancuso and Troy, 2020) at such high spatial resolution that cellular behaviour—in the form of bone remodelling sites—can be studied and the corresponding mechanical loading can be calculated. Subsequently, these methods can be used to investigate bone’s underlying mechanoregulated remodelling process, which may be the key to the development of patient-specific therapeutic or pharmacological interventions for various bone diseases.

Typically, when investigating bone mechanoregulation under controlled experimental conditions, micro-FE models disregard subject-specific variations in external loading conditions using simplified uniaxial compressive displacement boundary conditions (SC) (Schulte et al., 2013; Mancuso and Troy, 2020; né Betts et al., 2020; Malhotra et al., 2021). However, when investigating mechanoregulation in patients, variations in day-to-day external loading are more substantial due to habitual differences and patient-specific variability in the musculoskeletal system’s performance. Distinctive tensile forces and moments are applied to joints on a routine basis to stabilise under gravitational and other external loads and create unique loading patterns (Watkins, 2009). Consequently, to investigate mechanoregulation under day-to-day loading in a personalised medicine approach, patient-specific physiological loading patterns and boundary conditions need to be estimated (Galibarov et al., 2010; Yosibash et al., 2020).

In an effort to quantify in vivo loading patterns using biomechanical models, several load estimation algorithms have been developed. Artificial neural network-based approaches have been proposed (Garijo et al., 2014, 2017; Mouloodi et al., 2020) but lack interpretability, which is critical for moving to diagnostic use in patients to guide local therapeutic interventions. As a result, an algebraic method introduced by Christen et al. (2012) has been widely implemented to approximate the internal load history based on bone morphology (Christen et al., 2014; Badilatti et al., 2017; Synek et al., 2019; Cheong et al., 2020; né Betts et al., 2020). This algorithm superimposes and scales a finite number of loading cases until a target tissue load of homogeneous strains is found. Christen et al. (2012) demonstrated the capabilities of such a reverse-engineering approach using an extra-physiological tail-loading animal model, predicting the applied compressive loading in mouse caudal vertebra. However, the remaining signal inhomogeneity remained high, ranging between 20% and 67%, indicating that no homogeneous tissue load could be found (Christen et al., 2012). This suggests that only part of the bone structure may be load adapted. The actual in vivo load distribution might differ systematically from the homogeneous assumption in humans (Christen et al., 2016; Johnson and Troy, 2018) and mice (Christen et al., 2012). By modelling homogenised strain patterns, the conventional algorithm may reduce mechanical signal inhomogeneities that have been recognised as drivers for the mechanoregulated remodelling process in bone (Frost, 1987, 2003). Thus, this model’s assumptions may not be optimal and do not fully utilise all available information in time-lapsed data of longitudinal bone imaging studies.

This study had two goals. First, to derive an in silico-validated, robust, and specific method to estimate in vivo loading. Second, to apply this algorithm to examine in vivo mechanoregulation (Schulte et al., 2013) in humans and mice. We hypothesised that by extracting bone remodelling sites from time-lapsed imaging data, the relationship between bone formation in high-strain regions, quiescence in medium-strain regions, and resorption in low-strain regions could be used in a reverse-engineering optimisation approach to determine organ-level loads. We verified our mechanoregulated approach (MR) using three unique datasets and benchmarked it with an existing load history (LH) algorithm (Christen et al., 2012). First, to calculate sensitivity, specificity, and accuracy, MR and LH algorithms were applied to synthetic remodelling data derived from HR-pQCT images (Badilatti et al., 2016; Ohs et al., 2020a). Second, to test whether the algorithms are capable of predicting the loading conditions in a controlled experimental setup, both algorithms were applied to micro-CT scans of two groups of mice that had their caudal vertebra either loaded (8 N) or sham loaded (0 N) from a previous study (Scheuren et al., 2020b). Third, to assess the method’s fidelity in patients, MR and LH algorithms were applied to time-lapsed HR-pQCT scans and compared to patient-specific handgrip force measured using a dynamometer. Finally, to quantify the association between bone remodelling and mechanical stimulus, we derived a correct classification rate (CCR) (né Betts et al., 2020).

With the increasing prevalence of bone mechanoregulation studies, this work aimed to extend a previously developed load estimation algorithm (LH) (Christen et al., 2012) by allowing for tissue strain inhomogeneities in our mechanoregulated load estimation approach (MR). These localised differences in mechanical signal may drive bone’s remodelling response and help understand bone mechanoregulation. We provided validation for both algorithms using in silico-generated data, in vivo HR-pQCT images in humans, and micro-CT images in mice. These experiments indicate the portion of bone remodelling that can be attributed purely to mechanics and establish a baseline for futures studies evaluating mechanoregulation in patients.

Importantly, a combined in silico validation and in vivo verification, as shown in this study, has not yet been carried out. As such, algorithmic performance quantification was able to be carried out in human radius geometries and mice. Previous studies provided validation using in vivo mouse loading experiments (Christen et al., 2012). However, this did not enable the demonstration of algorithmic functionality for load directions other than uniaxial compression, such as those observed in the human distal radius and the mouse vertebra. The consistent results between our in silico and in vivo loading experiments indicate the validity of the MR algorithmic assumptions under diverse loading conditions. Corroborating the necessity for algorithmic validation in all six degrees of freedom, our in silico experiments identified possible performance deficits when applied to complex loading regimes. Despite using the inverse mechanoregulation rules of the advection simulation, MR’s in silico accuracy was not perfect for several reasons. First, only six selected time points (out of 40 simulated remodelling steps) were used to generate an in silico HR-pQCT scan series that reflected our in vivo data. Consequently, the inverse optimisation was challenged to recover loading from an iteratively adapted structure in a single step. Second, the advection simulation’s force-controlled setup caused slight differences in remodelling rates due to the initial anisotropy of the physiological load-adapted bone structure. Here, CS resulted in slightly higher average SED and BV/TV values by favouring bone formation compared to other load scenarios. In contrast to MR, the advection model limits the maximum bone formation rates, which may partially explain the performance deficits within this group.

Although MR’s performance was excellent for simulated adaptation, in vivo bone remodelling is not purely load-driven. Predicted in vivo loading patterns in the mouse model were in good agreement with a previous study (Christen et al., 2012). Compared to the dataset used by Christen et al. (2012), our LH results showed slightly larger moments while MR predictions were overall in good agreement with the previous study. Our LH results suggest a sizeable torsional component was induced in the caudal vertebra during daily activity, conflicting with the fact that the intervertebral discs limit the transmission of axial moments. Following the model proposed by Schulte et al. (2013), the intervertebral discs of the mouse FE analysis were modelled as stiff tissue (14.8 GPa), which may have resulted in slightly more uniaxial loading. The torsional moment may be affected by this assumption combined with the homogeneous strain simplification of LH as it was not detected using MR. LH-predicted in vivo forces in the distal radius model were consistent with a previous study (Christen et al., 2013); however, predicted moments varied by order of magnitude. Christen et al. (2013) used layers of soft-tissue padding at the proximal and distal ends, which may have resulted in further homogenisation of the strains throughout the radius. As such, this step may have limited the transmission of moment load at the interface between calcified tissue and soft tissue. When comparing our results with a cadaveric study investigating distal radius load during various wrist motions (Smith et al., 2018), we find similar load-to-moment ratios, indicating that additional padding may lead to an overestimation of momentum load. Finally, processes such as calcium homeostasis, wherein random bone remodelling may occur, will influence MR estimations. However, the findings of our mechanoregulation analysis reveal that the strain patterns overlap with the pattern observed by natural bone remodelling activity and can be used to estimate in vivo loading through our MR reverse engineering approach.

Our data also suggest that estimates in the distal radius may vastly vary from patient to patient. Despite the variance, an increase in loading was associated with increases in grip strength among patients. Such a relationship has been previously reported in cadaveric studies correlating grip strength with joint forces. In agreement with our results, models showed that approximately 26.3 N of force needs to be transmitted through the radius to obtain 10 N of grip strength (Putnam et al., 2000). Although this correlation was significant for loads estimated by MR and LH in the current study, this relationship was largely driven by single individuals with high grip strength. For future distal radius studies, grip strength should be considered as an inclusion parameter. Overall, our results indicate that the internal loads estimated by MR and LH are in good agreement with previous studies and can be linked to external factors such as grip strength in patients.

The principal algorithmic differences between approaches establish different future applications for LH and MR. MR prioritises remodelling sites, which are derived from two subsequent time-lapsed images. Accordingly, MR’s estimation is limited to the time frame between scans. LH estimates loading based on the bone morphology and is therby a cumulative estimate of all prior loadings (load history). Loading during immobilisation treatment (Lill et al., 2003; Clayton et al., 2009; Spanswick et al., 2021), exercise (Troy et al., 2020), or loading interventions (Hughes et al., 2018) may differ from a patient’s load history, which is defined by everyday and occupational activities. Thus, cumulative estimates of LH may be biased by the initial conditions. We showed that initial calibration of LH tends to improve differentiation between loading scenarios; however, this does not allow LH to achieve the same performance as MR. For the mouse loading experiment, this was evident in the delayed detection of significant differences between the loaded and control groups using LH compared to MR. For the present distal radius dataset, patients were skeletally mature adults and did not participate in a specific loading intervention. As a result, there was good agreement between MR and LH estimates. Note that the intact, contralateral radii used in the present study were taken from a patient cohort that had experienced a distal radius fracture. As such, loading in the unfractured arm may have increased, particularly in cases involving fracture of the dominant arm. The resulting change in day-to-day loading may explain slightly higher predictions of MR compared to LH throughout the study. While our results indicate that MR is more sensitive to changes in loading, the algorithm is also more affected by imaging bias than LH. By utilising two subsequent HR-pQCT images, MR is subject to higher noise levels, movement artefacts, and registration errors compared to LH (MacNeil and Boyd, 2008; Sode et al., 2011). LH may be more suited to mouse studies, which can assess lifetime changes, but not for the time frame of most clinical studies of antiresorptive therapies that often assess changes in BMD over a study duration of less than 2 years (Chen and Sambrook, 2012). Overall, our results have confirmed MR’s and LH’s capabilities for various applications using well-defined in silico loading and controlled experimental conditions. Accordingly, MR should be used when investigating designated time intervals in a longitudinal analysis and LH to assess the loading history in a cross-sectional fashion or when confronted with low image quality.

To quantify mechanoregulation, we have used a CCR similar to the approach described by né Betts et al. (2020). Here, we show that by using the boundary condition derived by MR, we achieve significantly higher CCR values than LH for simulated, physiological, and extra-physiological loading. Furthermore, our results indicate that these differences are more pronounced when an extra-physiological load was induced. Our results also show that using the simplified compressive boundary condition may be an acceptable choice when investigating trabecular bone mechanoregulation of the healthy human distal radius. However, Johnson and Troy (2018) have shown that this simplified compression boundary condition may alter cortical and trabecular loading sharing. Therefore, the authors caution that such a simplified boundary condition may not be adequate for future studies investigating cortical and trabecular bone mechanoregulation. Although our results indicate a higher trend in CCR for loads estimated by MR, we cannot entirely rule out the possibility that inherent parallels between mechanoregulation analysis and MR synthetically inflate CCR within human distal radius data. However, our analysis of an in vivo loading model has provided experimental ground truth showing that estimations by MR reflected experimental conditions properly in mice. Furthermore, our in silico validation showed that MR is highly sensitive, specific, and accurate. Overall, our results indicate that mechanoregulation tends to be higher when analysing physiological loading derived by MR and thrives on a wealth of extra-physiological loading. Interestingly, our results also show that simple compression is an adequate simplification for the in vivo loading environment in the distal radius considering current limitations. Furthermore, the results of our mechanoregulation analysis revealed a pronounced positive correlation between bone resorption and low strains for our mouse and a human model. This is in agreement with a previous study by né Betts et al. (2020) investigating mechanoregulation in a rodent femoral defect model, which indicated that mechanoregulated bone resorption mainly occurred within the distal and proximal fragments early during recovery. This relationship would indicate that osteoclastic activity may be more sensitive to local strain, and mechanoregulation may differ locally throughout the bone.

The proposed MR algorithm is subject to several limitations attributable to model assumptions as well as experimental and computational constraints. The performed multiclass ROC analysis weighted percent deviations in loading between forces and moments expressed as a reference force, equally. Where forces and moments may have a different impact on tissue level SED, the underlying in silico experiments (C, SX, SY, MZ, BX, BY, CS, CB, and 6DoF) were performed at equal loading magnitudes for each geometry, making this method a reasonable in silico performance measure. Regarding the animal experiments, the adjacent vertebra’s pinning procedure is limited in precision, and vibrations during the vertebra loading may create slight variations in loading direction and explain the observed higher variability in lateral bending. However, our results are comparable to a previous study (Christen et al., 2012) and represent the experimental setup sufficiently to provide validation for MR and LH. Regarding computational aspects, the method used to determine remodelling sites may include artefacts from scanning, such as beam hardening, motion artefacts, and partial volume effects or numerical inaccuracies of the image registration. However, in vivo micro-CT and HRp-QCT have been shown to have sufficient reproducibility for longitudinal bone structure assessment (Ellouz et al., 2014; Scheuren et al., 2020a). Additionally, MR used a Nelder–Mead optimiser that is not a true global optimisation algorithm and may converge in a local solution. However, in practise, it tends to work reasonably well for nonlinear, multimodal, inherently noisy functions. To further counteract this effect, we initialised the optimiser using a least-squares solution (as derived by LH). Future studies confronted with lower image quality may consider using Bayesian Global Optimisation techniques, which come at a higher computational cost but exhibit statistical methods, to address this problem. According to previous work, SED was used as a mechanical signal (Christen et al., 2012, 2013, 2014). More recent studies (né Betts et al., 2020; Malhotra et al., 2021) have identified an effective strain as a preferred candidate for bone mechanoregulation analysis using multi-density FE analysis. However, previous research has shown that these signals are strongly correlated (Pistoia et al., 2002; Ruimerman et al., 2005). Also, the FE model used was linear regarding material and geometry, and load cases were scaled and superimposed linearly during the optimisation to model the compounded loading effect. These simplifications would not capture any nonlinear behaviour or viscoelastic effects; however, only small linear-elastic deformations are expected to occur during day-to-day activity. Future studies may expand this model with increasing computational power and investigate nonlinear effects above yield strength that lead to bone failure (Schwiedrzik and Zysset, 2015).

We have shown that MR is an enhanced load estimation algorithm tailored for longitudinal bone remodelling studies, achieving high sensitivity, specificity, and accuracy in silico by employing acknowledged mechanoregulation principles. The combined in silico validation and in vivo verification approach presented in this study proved to be a powerful benchmarking tool for the development of time-lapsed bone imaging analysis methods. Moreover, our results indicate that future studies may use grip strength as a functional surrogate to verify estimated patient-specific physiological distal radius loads. Finally, our mechanoregulation analysis revealed considerable amounts of mechanically driven remodelling activity driven in human bone that may enable future studies to understand osteodegenerative disease.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

The studies involving human participants were reviewed and approved by Medical University of Innsbruck Ethics Committee. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by Swiss Federal Food Safety and Veterinary Office.

MW, NO, MB, RM, and CC: study design. MW and FM: study conduct. MW and MB: data collection. MW, NO, and FM: data analysis. MW, FM, NO, RM, and CC: data interpretation, revising manuscript content, and approving final version of manuscript. MW and CC: drafting manuscript. CC: takes responsibility for the integrity of the data analysis. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.