Cortical Thickness Adaptive Response to Mechanical Loading Depends on Periosteal Position and Varies Linearly With Loading Magnitude

Introduction

Osteoporosis is a degenerative bone disease characterised by long-term bone loss and fragility (Black and Rosen, 2016). To counter osteoporosis, several drugs have been developed to either reduce or reverse the bone loss process. Despite the success in reducing the degeneration of osteoporosis, drug treatments can have significant side effects, and the positive effect on the bone mass is often lost upon discontinuation of the drug dosing regimen (McClung, 2016; Minisola et al., 2019). Pharmacologically, drug treatments such as PTH are generalised therapies and do not target specific bones. Exercise, on the other hand, has been identified as a safe alternative to restore bone mass (Bliuc et al., 2013; Ebeling et al., 2013; Beck et al., 2017); mechanical loading interventions can act as a potent anabolic stimulus with the ability to strategically restore bone mass in regions of bone that undergo significant loading, both in animal models and humans (Ozcivici et al., 2010).

Bone tissue adapts its mass and structure to the habitual mechanical loading environment (Rubin and Lanyon, 1985; Ozcivici et al., 2010; Pivonka et al., 2018). Several animal loading models have been developed to investigate the relations between the applied mechanical load, the changes in bone mass, and the bone cells involved in mechano-transduction (Meakin et al., 2014; Javaheri et al., 2019). Among these models, the mouse tibia loading model is commonly used to assess both trabecular and cortical bone adaptation responses (De Souza et al., 2005; Sugiyama et al., 2010, 2012). In this model, the tibia is subjected to cyclic, compressive load, while the contralateral tibia serves as an internal control. Common metrics used to assess bone adaptation to mechanical loading consider global morphological variations at either the entire bone level (e.g., bone volume change), or at the entire bone slice level (e.g., cross-sectional area change, moment of area change). A comprehensive study by Sugiyama et al., 2012 (Sugiyama et al., 2012) explored the influence of peak dynamic load on bone adaptation. The primary focus was to analyse the effect of peak dynamic loads (ranging from 0 to 14 N) on changes in cortical area (ΔCt.Ar), determined through μCT endpoint imaging. They concluded that changes in cortical bone cross-sectional area are linearly related to the peak applied load. However, bone adaption is a local (i.e., site-specific) phenomenon governed by the local strain (Fritton et al., 2005; Razi et al., 2015). As such, metrics operating on the entire bone or on a slice level are unable to provide detailed insights into a load-adaptation response law. While the study of Sugiyama et al. did observe site-specific adaptation, the load-adaptation response was not explored this in a quantitative manner. A more detailed evaluation of bone’s adaptive response to local strain can be obtained by analysing the local cortical thickness change (ΔCt.Th).

To this end, several studies have explored the local cortical thickness variation (ΔCt.Th) (Halloran et al., 2002; Stadelmann et al., 2011; Sugiyama et al., 2012; Galea et al., 2015; Birkhold et al., 2016; Roberts et al., 2020), commonly using a minimum distance metric (i.e., the shortest distance between periosteal and endosteal surfaces) (Hildebrand and Rüegsegger, 1997; Bouxsein et al., 2010). Pereira et al. (2015) used the same method to analyse ΔCt.Th but instead considered spatially discrete locations, reporting ΔCt.Th in a polar coordinate system around the centroid. This technique provided promising results for the majority of the tibial cross-sections analysed; however, it is inadequate for bony protrusions such as the tibial ridge. Furthermore, while new tissue forms normal to the bone surface (Graham et al., 2012; Pereira et al., 2015; Zhang et al., 2019), the use of a minimum distance technique does not accurately capture ΔCt.Th when the direction of adaptation (i.e., normal to the surface) is highly offset from the radial direction. Similarly, radial coordinates create issues when determining periosteal and endosteal edges, where in some cases up to four cortical intersection points can be identified for a given radial direction (Bab et al., 2007).

To account for the irregular shape of the mouse tibia, this study proposes a new technique for measuring cortical thickness variations. A novel image post-processing algorithm was developed to allow the calculation of the local ΔCt.Th around the perimeter of the tibia using a combined minimum distance and normal distance approach. The experimental results of Sugiyama et al. (2012) were re-analysed to quantify local cortical thickness changes and their association to the peak load applied. The analysis was conducted for four commonly studied cross-sections in the mouse tibia loading model (i.e., 25, 37, 50, 75%). Furthermore, mechanical analysis using beam theory was performed in order to relate the obtained cortical thickness changes to the local mechanical loading environment and identify regions of high and low strains, respectively.

Materials and Methods

The endpoint imaging data used in this study was previously reported by Sugiyama et al. (2012). As such, we have provided a brief summary of the experimental design and imaging process here; for a more complete description see (Sugiyama et al., 2012). Following this, a detailed description of the newly developed image processing algorithm used to extract local Δ Ct.Th measurements of tibial cross sections at selected regions is presented.

Experimental Design

A total of 48 female C57BL/6 mice were divided evenly into eight groups, with each group assigned to one of eight peak load magnitudes (F = 0, 2, 4, 6, 8, 10, 12 or 14 N) (Sugiyama et al., 2012). For the purpose of our study, the F = 14 N loading case was excluded due to the formation of woven bone in several animals. Each mouse was subjected to a right sciatic neurectomy at 17 weeks of age, in order to minimise the natural loading in their right tibiae (i.e., muscle contraction forces) and simulate a condition of mechanical disuse. From day 5 after neurectomy, every second day, and for two weeks, the right tibia of each mouse was subjected to external mechanical loading. A non-invasive servo hydraulic loading machine applied 40 cycles of intermittent loading, with each cycle consisting of: (i) 0.5 N static preload, (ii) 500 N/s ramp up to target peak load, (iii) a 0.05 second hold at peak load, (iv) −500 N/s ramp down to static preload, (v) 10 s rest interval. This has been shown to significantly stimulate loading-related bone gain (Rubin and Lanyon, 1984; Fritton and Rubin, 2001; Robling et al., 2001; Srinivasan et al., 2002; De Souza et al., 2005; Sugiyama et al., 2010, 2011; Moustafa et al., 2012). The left tibia of each mouse was used as contralateral control (Sugiyama et al., 2010; McKenzie and Silva, 2011). At day 21 after neurectomy, the mice were sacrificed, and both left and right tibiae were scanned using μCT imaging. Whole tibiae were imaged using the SkyScan 1172 (SkyScan, Kontich, Belgium) with an isotropic resolution of 4.78 μm. An X-ray voltage of 50 kV was applied, with 0.5 mm aluminium filtration. The scans were over 180 degrees with a 0.5-degree rotation step.

Beam Theory Analysis of Tibia

The mechanical analysis presented in this paper aims to link the strains in the cortical cross section to the observed thickness changes. We assume that the tibia represents a slender beam structure and, consequently, can be analysed using Euler-Bernoulli beam theory (Hjelmstad, 2005; Bauchau and Craig, 2009; Buenzli et al., 2013; Lerebours et al., 2016; Trichilo, 2018; Ashrafi et al., 2020). The purpose of this analysis is not to provide a direct link between strain magnitudes and the adaptive response, but rather to identify compressive and tensile regions of strain and to observe general trends of strain magnitude across a given cross-section.

The load F was assumed to act on the tibial plateau in the z-direction between the tibial condyles; this location was previously suggested from strain gauge studies (Pickering et al., 2021). In a particular cross section (z), F induces a normal force (F = N) and bending moments M_x ( = F⋅I_y) and M_y ( = F⋅I_x), where I_y and I_x represents the distance of the load F to the x and y axis respectively. Knowing the internal beam quantities once can calculate the axial strain according to:

$\begin{array}{c}\mathrm{\epsilon }(x,y)=\frac{1}{E}\mathrm{\sigma }(x,y)=\frac{1}{E}\left(\frac{N}{A}-\frac{M_y{I}_x+M_x{I}_{xy}}{I_x{I}_y-I_{xy}^2}x+\frac{M_x{I}_y+M_y{I}_{xy}}{I_x{I}_y-I_{xy}^2}y\right)\end{array}(1)$

where I_x and I_y are the second moments of area with respect to the x- and y-axis, respectively and I_xy is the product moment of area. In Eq. 1 bone was assumed to be a linear-elastic material with a Young’s modulus E = 17 GPa (Kohles et al., 1997). A maximum load of F = 12 N was applied to aid in differentiation of strain magnitudes, aligning with the peak load used in the experimental protocol.

To compute the second moments of area I_x, I_y, and I_xy of each cross-section, a customised algorithm to automatically segment μCT images was developed in MATLAB. μCT images were first binarised using Otsu’s thresholding method. A filter was then applied to the images to close any small holes in the tibial cross-section (representing blood vessel channels), and to smooth the tibial boundaries. The second moments of area were calculated using parallel axis theorem, treating each white pixel (i.e., bone) in the image as a square of area 22.84 μm².

Image Post-processing Algorithm

Each stack of tibial μCT images was normalised along the proximal-distal direction of the tibia (i.e., z-axis), with z = 0% referring to the most proximal slice and z = 100% referring to the most distal slice. In this study, the response to mechanical loading was analysed on a single cross-sectional slice taken from the z = 25% (proximal), z = 37% (proximal-middle), z = 50% (middle), and z = 75% (distal) locations of the tibia. Note that we also performed the following methodology over a representative stack of images spanning approx. 0.5 mm of the tibia (±0.25 mm from selected slice), as has been commonly done in previous works (see Sugiyama et al., 2012). A comparison of the results from using a single slice and the representative stack can be found in Supplementary Figure 1; measurements in a single slice were found to not differ significantly from the representative stack. The selected μCT images were grouped based on peak load applied, cross-section analysed and control/loaded tibia.

Following the binarisation process described above, pixels along the periosteal and endosteal envelopes were identified and mapped into an array. In order to compare the thickness along the periosteum between different limbs, periosteal position (P_per) distributions were aligned across all tibiae at a given z cross-section through the location of a characteristic point (i.e., pixel) on the tibial periosteum (P_per = 0). This characteristic point was identified as the intersection between the tibia periosteum and the line connecting the tibia and fibula centroids (Figure 1A). Starting at P_per = 0 and following a clockwise direction, the pixels along the periosteal surface were re-arranged and normalised between 0 and 1. In the case of the distal cross-section (z = 75%) where the fibula is absent, a faux fibula centroid was projected onto the plane from the 50% section.

FIGURE 1

Figure 1. Measurement methodology used to analyse the cortical thickness of mice tibia (z = 25% section). (A) Pre-processing – the cortical boundaries are located (periosteum and endosteum), and a characteristic periosteal point is determined as the intersection of a line connecting the centroids of the tibia and fibula (red dot). The periosteal length (perimeter) is normalised between 0 and 1. Ct.Th measurements at each periosteal point are taken clockwise around the tibial perimeter. (B) Thickness measurement case 1: minimum distance method (red) and surface-normal blue) measurements around the cortex. Similar results with minimum distance providing shorter measurements. (C) Thickness measurement case 2: thickness measurements along the tibial ridge, showing major differences between measurement methods with tangent-normal providing shorter measurements. (D) Combination of the two measurement methods, selecting the smallest distance determined by either measurement method, to create the most representative cortical thickness distribution of the tibial image. Combined Results are then filter using a 2nd order Butterworth filter.

When defining the local cortical thickness two types of measurements were used, as shown in Figures 1B,C: (i) a minimum distance measurement (shown in blue) and (ii) a perpendicular distance measurement (shown in red). The minimum distance method measured the distance between each periosteal pixel and the nearest endosteal pixel. The normal method measured the distance from each periosteal pixel to the next cortical edge (endosteal or periosteal) along a line perpendicular to the periosteum. In cross-sectional regions with approximately constant curvature of periosteum and endosteum, both techniques provided a similar result, as shown in Figure 1B. The more distally located cross-sections (z = 50 and 75%) confirm this trend for the majority of the periosteal surface. However, in regions with large curvature changes such as the tibial ridge (i.e., z = 25 and 37%), a large discrepancy between the two measurement techniques was observed, as highlighted in Figure 1C. In order to generate a thickness measurement which best represented bone adaptation in these sections [i.e., normal to the surface (Graham et al., 2012)], both thickness measures were calculated for each periosteal pixel and the smaller of the two measurement values was used to define the representative local cortical thickness Ct.Th (Figure 1D). This result was then filtered using a 2nd order lowpass Butterworth filter to remove high-frequency noise due to the measurement combination technique.

Since the results across images were different in length due to the variability between animals and the adaptation process of the loaded limbs, each Ct.Th distribution was re-sampled so that P_per contained n = 750 periosteal points. A further consideration was made when comparing cortical thickness distributions of the loaded and control limbs. As mentioned, the loaded limb presented a longer periosteum due to the adaptation process, therefore, to ensure an accurate comparison of the same cortical regions, a further alignment step was required. For the approximately circular cross-sections (z = 50 and 75%), cross-covariance was used to circularly shift one of the two thickness signals, to maximise the alignment with the other one. In the cases of the z = 25 and 37% cross-sections, where the growth/resorption along the tibial ridge had a significant effect on the alignment of P_per points between loaded and control tibia, a customised re-sampling methodology was developed. Four common peaks and/or troughs in the Ct.Th measurement distribution were identified in all the limbs analysed that correlated with key bony features, e.g., the tip of the tibial ridge. Thickness measurements between these key-points were re-sampled based on a fixed number of points, resulting in an optimal and consistent alignment between the peaks for all slices at that particular cross-section. This process was repeated for each mouse limb (right and left tibiae). It should be noted that the thickness along the periosteum was measured for both right and left limbs starting from P_per = 0 and following a clockwise direction in the cross-section. To be able to make left vs. right comparisons, all the left limb signals needed to be reversed. The thickness distributions were then compared between the right and the left tibiae of a mouse at each cross-section (z), for all considered loading conditions. The relative change in cortical thickness (ΔCt.Th) for each periosteal point was calculated as:

$\text{Δ}Ct.Th(P_{per})[\%]=\frac{Ct.Th{(P_{per})}_{right}-Ct.Th{(P_{per})}_{left}}{Ct.Th{(P_{per})}_{left}}\cdot 100(2)$

where P_per identifies the periosteal position at which the cortical thickness is evaluated at the right and left tibia. Eq. 2 is an extension of the equation used for calculating cortical area changes (see Sugiyama et al., 2012 for details) with respect to considering localised cortical thickness changes.

Statistical Analysis

Mean values ($\overline{\mathrm{\Delta }\mathrm{Ct}.\text{Th}}$) and standard deviations (SD) of ΔCt.Th were calculated across the six specimens within a loading condition at each point P_per. For simplicity of notation, the symbol ΔCt.Th will denote mean cortical thickness changes throughout the rest of the manuscript. The results were evaluated through Student’s paired t-tests for each loading condition per tibial cross section, considering the link between the local thickness change at each periosteal surface position and the mechanical loading. For this investigation, a p-value <0.05 was considered statistically significant.

Results

For clarity, results will first be presented for the middle region (z = 50%) as this section consists of an approximately circular cross-section and embodies the greatest strains. Results for the remaining three sections will be presented thereafter.

Figure 2 shows the results of the Student’s paired t-tests for the 12 N loading case in the middle cross-section of the tibia. From Figure 2A one can see that the mean ΔCt.Th reaches values greater than 50% at selected positions along the periosteum (P_per ≈ 0.02), whereas other regions show zero or negative thickness change (0.3 < P_per < 0.35 and 0.7 < P_per < 0.80). Furthermore, mice within each loading group showed variable response to mechanical loading, with the standard deviation being approximately ±10% at the majority of periosteal surface. The results of the paired t-test indicate that local changes in cortical thickness were statistically significant (p < 0.05) over a large portion of the periosteal surface (Figure 2B), and the only regions demonstrating no significance were those with near zero cortical thickness change.

FIGURE 2

Figure 2. Local cortical thickness variation for the 12 N loading case at the middle tibial region (z = 50%): (A) Mean (black line) and standard deviation (Shaded area) along the periosteum ΔCt.Th ± SD vs P_per and (B) p-value using a student-paired t-test at each periosteal position (statistical significance indicated by values below the dashed line, i.e., p < 0.05).

Figure 3 shows the adaptation response along the periosteal surface in the z = 50% section for all peak loads investigated. For low peak loads (i.e., 4 N and lower), resorption was observed at most locations around the periosteal surface (negative ΔCt.Th, Figure 3A). At peak loads of 6 N and higher, positive ΔCt.Th (i.e., bone formation) begins to show along the periosteum. The region on the periosteum with bone formation increases with the increasing peak load. Note that the ΔCt.Th vs. F response is quasi linear, while the slope of the response depends on the periosteal position.

FIGURE 3

Figure 3. (A) Mean cortical thickness changes (ΔCt.Th) along the normalised periosteal position at the middle tibial cross-section (z = 50%) for all loading cases. Regions of interest have been identified at P_per = 0.05 (posterior surface), 0.33 (lateral), 0.50 (anterior), and 0.75 (medial). (B) Mean cortical thickness changes at selected regions of interest across all loading conditions. (C) Finite element results (F = 10 N) within the 50% cross-section.

At approximately P_per = 0.33 and P_per = 0.75, consistent resorption was observed with little to no dependence on the applied load. In contrast, periosteal locations at approximately P_per = 0.05 and P_per = 0.5 show a large dependence on load magnitude, exhibiting near maximum changes in ΔCt.Th. To explore the load dependency further, ΔCt.Th has been reported in Figure 3B as a function of applied load for the four periosteal locations identified above, i.e. P_per = 0.05, 0.33, 0.50, 0.75. At locations P_per = 0.33 and P_per = 0.75 a nearly constant reduction of ΔCt.Th, independent of the magnitude of the applied peak load, is noted. At locations P_per = 0.05 and P_per = 0.50, a quasi-linear relationship between load magnitude and ΔCt.Th is observed. The P_per = 0.05 location (i.e., posterior-lateral surface) was noted to experience greater ΔCt.Th when compared to P_per = 0.50 (i.e., anterior-medial surface).

To test if there is a correlation between the thickness change along the periosteum and the axial strain (ε_z) encountered in the cross section, beam theory was used to calculate ε_z in the middle cross-section (Figure 3C). Comparing the ΔCt.Th distribution (Figure 3A) and the axial strain (Figure 3C), apparent trends of the load-adaptation response can be observed; to investigate this further, four locations around the periosteum [posterior (P_per = 0.05), lateral (P_per = 0.33), anterior (P_per = 0.5), medial (P_per = 0.75)] were extracted and explored in Figure 3B. The posterior region experienced a higher strain magnitude compared to the anterior region (−8,198 and 6,075 με, respectively), coinciding with a higher ΔCt.Th along the periosteum in the same region. Likewise, the load independent locations on the periosteum (P_per = 0.33 and P_per = 0.75) appear to coincide with near zero axial strain. This region is commonly referred to as the neutral bending axis or neutral axis.

Figure 4 shows the mean cortical thickness changes, i.e., ΔCt.Th curves of the proximal (z = 25%), proximal-middle (z = 37%), and distal (z = 75%) cross-sections of the tibia, as well as the respective strain distribution ε_z generated by a 12 N load for the corresponding cross sections. The proximal and proximal-middle cross-sections show similar trends in ΔCt.Th to the middle region discussed above (Figures 4A,B), with bone gain or loss responses occurring at different periosteal positions. In the proximal cross-section, the maximum bone formation response occurs in the tensile region of the cross-section (0.35 < P_per < 0.65). Cortical growth response to compressive loading produced up to a 55% increase in cortical thickness at P_per = 0.9 for 12 N peak load. The proximal-middle cross-section follows similar trends to the proximal one. Bone gain is higher in the tensile region (maximum change of 80% at P_per = 0.49 for 10 N peak load) than it is in compression (maximum 53% at P_per = 0.96 for 12 N peak load).

FIGURE 4

Figure 4. (A–C) Mean cortical thickness changes (ΔCt.Th) vs. normalised periosteal position (P_per) at the: (A) proximal, (B) proximal-middle, and (C) distal tibial cross-section for all different loading cases. Vertical dashed lines represent the approximate P_per position of the expected neutral axis (NA). (D–F) Beam theory results for the: (D) proximal, (E) proximal-middle, and (F) distal tibial cross-sections. Dashed lines represent the approximate physical locations of the neutral axis (ε = 0).

In both the proximal and proximal-middle cross-sections, bone resorption was observed around the medial neutral axis under all loading conditions (0.25 < P_per < 0.33 in the proximal region, 0.15 < P_per < 0.3 in the proximal-middle region). The lateral side of the neutral axis was observed to show smaller rates of resorption with increased load (0.6 < P_per < 0.7 in the proximal region, 0.65 < P_per < 0.75 in the proximal-middle region). Under the 12 N load, the proximal-middle region shows no loss of bone at the lateral side of the neutral axis.

The adaptive response in the distal cross-section of the tibia differs significantly from the other three regions (Figure 4C). Bone gains are the lowest of all the four investigated cross-sections, showing a maximum cortical thickness increase of 21% at P_per = 0.17. The anterior surface (0.35 < P_per < 0.65) shows that the cortical surface remains relatively stable with zero change to cortical surface for F ≥ 4 N. Significant bone gain was observed in the posteromedial section of the cross-section (0 ≤ P_per < 0.35) for F ≥ 6 N, whereas the posterolateral section experienced resorption across all loading conditions. The neutral axis of the distal region cannot be clearly identified from the obtained thickness results.

The results of the beam analysis revealed peak tensile strains of 2,583, 4,187, and 964 με at P_per ≈ 0.5 for the proximal (Figure 4D), proximal-middle (Figure 4E), and distal (Figure 4F) cross-sections, respectively. Peak compressive strains were −6,875, −9,994, and −3,381 με at P_per ≈ 0.05. This is consistent with the middle region, where higher strain magnitudes were observed on the posterior side when compared to the anterior. It should also be noted that the majority of the distal cross section is under compressive loading with only a small tensile region (Figure 4F).

Discussion

Bone adaptation is a local phenomenon. This has been demonstrated previously along the axis of the tibia or in discrete segments of the tibia (Sugiyama et al., 2010; Lu et al., 2016; Galea et al., 2020; Roberts et al., 2020; Scheuren et al., 2020). However, this work is the first to demonstrate the link between magnitude of axial load and adaptive response around the periosteum of the tibia at a cross-section level. Our results clearly show that the cortical thickness change around the periosteal surface varies linearly with loading magnitude. Furthermore, the slope of this adaptive response depends on the periosteal position.

Quantifying the adaptive response through ΔCt.Th provides important insights into understanding the localised changes compared to other metrics (e.g., ΔCt.Ar). Shown in Figure 3B, the bone’s adaptive response was found to have a quasi-linear relationship between load magnitude and bone formation. This observation supports previous findings, such as those of Sugiyama et al. (2012) who found ΔCt.Ar increased linearly with the applied load. However, the use of such metrics (i.e., non-localised) mask the true magnitude of adaptation. Under a peak load of 12 N at the 50% section, Sugiyama et al. (2012) recorded an average ΔCt.Ar of 15.5 ± 2.1% (Sugiyama et al., 2012). However, for the same load, at the same section, we show that ΔCt.Th can vary between −20 and +60%; this suggests that broad metrics such as ΔCt.Ar are insufficient to fully describe the adaptation response to mechanical stimuli.

Many studies have shown a link between local adaptive response and local strain magnitude (Fyhrie and Carter, 1986; Robling et al., 2006; Webster et al., 2012; Schulte et al., 2013; Lambers et al., 2015; Pereira et al., 2015; Carriero et al., 2018; Tiwari et al., 2018). However, none of these studies have looked at the effect of different loading magnitudes on the observed local changes. Here we compared the local adaptive response as a function of loading magnitude in different cortical cross-sections. Focusing on the 50% cross-section, the posterior surface (P_per = 0.05) experienced larger bone formation as a function of peak load (i.e., larger slope of the F vs. ΔCt.Th curve) compared to the anterior surface (P_per = 0.5), as shown in Figures 3A,B. The increased response in this region correlates to the strain magnitude resulting from combined bending and axial load. The peak compressive strain is greater than the peak tensile strain. As such, the compressive surface experiences a larger bone formation response. The same trend is seen in the 25 and 37% cross sections, shown in Figures 4A,B.

Regions of low strain (i.e., those near the neutral axis) experienced bone resorption. In the 50% cross-section, the medial and lateral surfaces (P_per ≈ 0.75, 0.33, respectively) experienced resorption, independently of the load applied. This finding stands in contrast to the study performed by Pereira et al. (2015). Most notably, in non-neurectomised C57BL/6 mice, adaptation to loading of F = 13 N was all positive, i.e., no bone resorption occurred (Pereira et al., 2015; Trichilo, 2018). One explanation for the observed bone loss in the neurectomised mice is the fact that muscle action comprises a significant portion of the habitual strain state. Without the continual influence of muscle activation to maintain mechanical homeostasis at the neutral axis, these regions will undergo bone resorption to readjust to their new habitual state.

This same trend of resorption around the neutral axis is observed in the 25 and 37% sections. In the case of the 75% cross-section, however, a more general trend of resorption is observed, with up to half of the cross-section experiencing resorption under the highest peak loads. Strains in this cross-section were noticeably smaller than the strains in the other three. This might be related to this cross-section of the tibia being aligned with the longitudinal axis, resulting in low bending moments, and thus axial compression dominates.

The bone adaptive response observed in Figures 3A, 4A,B can be further viewed through the lens of beam theory and second moment of area. A higher second moment of area leads to lower overall strain; the most efficient way to increase the second moment of area is by adding new material furthest away from the neutral axis. In doing so, bone maximises its strength while optimising the distribution of mass. This is a demonstration of Wolff’s law (Wolff, 1893). While a detailed mechanical analysis of the tibia using finite element analysis is beyond the scope of this paper, we are confident that the utilised beam theory predicted the location of the null axis and peak strains in the cross section well. Estimation of the exact magnitude of peak strain as a function of load is not the focus of the current paper.

Early understandings of bone’s adaptive response suggested a range of strain levels which would not elicit an adaptive response, often referred to as the lazy zone (Carter, 1984; Huiskes et al., 1987). More recent studies have suggested that this region is non-existent in both animal models and in human tissue (Sugiyama et al., 2012; Ellman et al., 2013; Schulte et al., 2013; Christen et al., 2014). Likewise, in this work, a lazy zone is not observed. In line with previous studies, the results presented here suggest that no lazy zone exists in adaptive bone (re)modelling.

A limitation of this study is that the measurement technique does not determine if adaptation has occurred on the periosteal or endosteal surface, rather it only determines the net thickness change. Furthermore, the neurectomy performed on the mice removes habitual loading and induces some amount of resorption, with loading inducing additional bone formation subsequently. Due to this phenomenon, we are unable to quantitatively determine the amount of new bone material formed or resorbed on each surface, only the total difference after completion of the experiment. Longitudinal imaging would provide a significant benefit in this regard. Comparing a single limb at different time points, differences on both the endosteum and periosteum could be tracked to provide deeper insights into the mechanisms of bone adaptation.

Effects of loading on neurectomised vs. non-neurectomised mice were also not considered in the original study. While the left limb was left intact, performing a sciatic neurectomy on the right limb may have affected the gait of the mouse, potentially altering the habitual strains experienced in the healthy limb. Without such a control, it is difficult to answer several questions such as how bone adapts to mechanical loads from a standard habitual state (i.e., healthy gait), what are the bone loss effects of neurectomy, and what are the differences in adaptation response between healthy versus mechanically deficient (i.e., neurectomised) mice. Answering these questions would help to provide a more complete understanding of bone adaptation and should be investigated in future studies.

In this paper, we presented a novel image processing algorithm to measure cortical thickness of the mouse tibia loading model and compared the results across several loading magnitudes. We identified that discrete locations around the periosteum were shown to follow a quasi-linear cortical thickness adaptation response with increased loading, while points at areas of near-zero strain (i.e., neutral bending axis) experienced resorption regardless of loading magnitude; the correlation between strain and bone formation was shown to follow the adaptation principles of Wolff’s Law.

The ultimate purpose of animal adaptation studies is to derive mechanistic insights into the link between applied mechanical loads and the observed organ- or tissue-scale changes of (cortical) bone. The work conducted here has established a statistically significant association of mechanical loading and bone adaptation responses in discrete periosteal regions of cortical bone. The fact that these regions also experienced high compressive and tensile strains obtained from beam theory provides confidence that a mechanistic relationship exists between a particular mechanical quantity (such as principal strain, strain energy density, etc.) and the local cortical thickness changes. These findings may be useful in the development of treatments that aim to increase bone strength, informing specific mechanical loading routines that would provide targeted bone formation in areas of high fracture risk. Results we have obtained here will help develop novel bone adaption algorithms which are able to predict cortical thickness changes which is the scope of a future study.

Data Availability Statement

The data analyzed in this study is subject to the following licenses/restrictions: Data used in this study is owned by the University of Bristol. Requests to access these datasets should be directed to PP, peter.pivonka@qut.edu.au.

Ethics Statement

The animal study was reviewed and approved by University of Bristol.

Author Contributions

CM: conceptualisation, development of algorithm to measure cortical thickness, beam theory analysis, collection of results, formal analysis, methodology, writing-original draft, and writing, review, and editing. ST: development of algorithm to measure cortical thickness, collection of preliminary results, and review and editing. EP: conceptualisation, beam theory analysis, data interpretation, and review and editing. SM: beam theory analysis and review and editing. PD and LM: data curation and review and editing. PP: conceptualisation, data curation, data analysis, funding acquisition, investigation, methodology, project administration, supervision, writing-original draft, and writing, review, and editing. All authors contributed to the article and approved the submitted version.

Funding

CM acknowledged the support of an APA scholarship by Queensland University of Technology. ST acknowledged the support of an international Ph.D. scholarship from the University of Melbourne. LM was supported by a Veterinary Training Fellowship provided by the Wellcome Trust (092045/Z/10/Z). PP and SM acknowledged the support from the Australian Research Council (IC190100020; DP180103146; FT180100338).

Conflict of Interest

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2021.671606/full#supplementary-material

References

Ashrafi, M., Gubaua, J. E., Pereira, J. T., Gahlichi, F., and Doblaré, M. (2020). A mechano-chemo-biological model for bone remodeling with a new mechano-chemo-transduction approach. Biomech. Model. Mechanobiol. 19, 2499–2523. doi: 10.1007/s10237-020-01353-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Bab, I., Hajbi-Yonissi, C., Gabet, Y., and Müller, R. (2007). Micro-Tomographic Atlas of the Mouse Skeleton. (Boston, MA: Springer), 171–181. doi: 10.1007/978-0-387-39258-5_16

CrossRef Full Text | Google Scholar

Bauchau, O. A., and Craig, J. I. (2009). Structural Analysis. Dordrecht: Springer Netherlands, doi: 10.1007/978-90-481-2516-6

CrossRef Full Text | Google Scholar

Beck, B. R., Daly, R. M., Singh, M. A. F., and Taaffe, D. R. (2017). Exercise and sports science Australia (ESSA) position statement on exercise prescription for the prevention and management of osteoporosis. J. Sci. Med. Sport 20, 438–445. doi: 10.1016/j.jsams.2016.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Birkhold, A. I., Razi, H., Duda, G. N., Weinkamer, R., Checa, S., and Willie, B. M. (2016). The periosteal bone surface is less mechano-responsive than the endocortical. Sci. Rep. 6:23480. doi: 10.1038/srep23480

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, D. M., and Rosen, C. J. (2016). Clinical practice. Postmenopausal osteoporosis. N. Engl. J. Med. 374, 254–262. doi: 10.1056/NEJMcp1513724

PubMed Abstract | CrossRef Full Text | Google Scholar

Bliuc, D., Nguyen, N. D., Nguyen, T. V., Eisman, J. A., and Center, J. R. (2013). Compound risk of high mortality following osteoporotic fracture and refracture in elderly women and men. J. Bone Miner. Res. 28, 2317–2324. doi: 10.1002/jbmr.1968

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouxsein, M. L., Boyd, S. K., Christiansen, B. A., Guldberg, R. E., Jepsen, K. J., and Müller, R. (2010). Guidelines for assessment of bone microstructure in rodents using micro-computed tomography. J. Bone Miner. Res. 25, 1468–1486. doi: 10.1002/jbmr.141

PubMed Abstract | CrossRef Full Text | Google Scholar

Buenzli, P. R., Thomas, C. D. L., Clement, J. G., and Pivonka, P. (2013). Endocortical bone loss in osteoporosis: the role of bone surface availability. Int. J. Numer. Method Biomed. Eng. 29, 1307–1322. doi: 10.1002/cnm.2567

PubMed Abstract | CrossRef Full Text | Google Scholar

Carriero, A., Pereira, A. F., Wilson, A. J., Castagno, S., Javaheri, B., Pitsillides, A. A., et al. (2018). Spatial relationship between bone formation and mechanical stimulus within cortical bone: combining 3D fluorochrome mapping and poroelastic finite element modelling. Bone Rep. 8, 72–80. doi: 10.1016/j.bonr.2018.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Carter, D. R. (1984). Mechanical loading histories and cortical bone remodeling. Calcif. Tissue Int. 36(Suppl. 1), S19–S24. doi: 10.1007/BF02406129

PubMed Abstract | CrossRef Full Text | Google Scholar

Christen, P., Ito, K., Ellouz, R., Boutroy, S., Sornay-Rendu, E., Chapurlat, R. D., et al. (2014). Bone remodelling in humans is load-driven but not lazy. Nat. Commun. 5:4855. doi: 10.1038/ncomms5855

PubMed Abstract | CrossRef Full Text | Google Scholar

De Souza, R. L., Matsuura, M., Eckstein, F., Rawlinson, S. C. F., Lanyon, L. E., and Pitsillides, A. A. (2005). Non-invasive axial loading of mouse tibiae increases cortical bone formation and modifies trabecular organization: a new model to study cortical and cancellous compartments in a single loaded element. Bone 37, 810–818. doi: 10.1016/j.bone.2005.07.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebeling, P. R., Daly, R. M., Kerr, D. A., and Kimlin, M. G. (2013). An evidence-informed strategy to prevent osteoporosis in Australia an outline of the building healthy bones throughout life white paper. Med. J. Aust. 198, 90–91. doi: 10.5694/mja12.11363

CrossRef Full Text | Google Scholar

Ellman, R., Spatz, J., Cloutier, A., Palme, R., Christiansen, B. A., and Bouxsein, M. L. (2013). Partial reductions in mechanical loading yield proportional changes in bone density, bone architecture, and muscle mass. J. Bone Miner. Res. 28, 875–885. doi: 10.1002/jbmr.1814

PubMed Abstract | CrossRef Full Text | Google Scholar

Fritton, J. C., Myers, E. R., Wright, T. M., and Van Der Meulen, M. C. H. (2005). Loading induces site-specific increases in mineral content assessed by microcomputed tomography of the mouse tibia. Bone 36, 1030–1038. doi: 10.1016/j.bone.2005.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Fritton, S. P., and Rubin, C. T. (2001). “In vivo measurement of bone deformations using strain gauges,” in Bone Mech. Handbook, 2nd Edn, ed. S. C. Cowin (Milton Park: Taylor & Francis), 8-1-8–42. doi: 10.1201/b14263-11

CrossRef Full Text | Google Scholar

Fyhrie, D. P., and Carter, D. R. (1986). A unifying principle relating stress to trabecular bone morphology. J. Orthop. Res. 4, 304–317. doi: 10.1002/jor.1100040307

PubMed Abstract | CrossRef Full Text | Google Scholar

Galea, G. L., Delisser, P. J., Meakin, L., Price, J. S., and Windahl, S. H. (2020). Bone gain following loading is site-specifically enhanced by prior and concurrent disuse in aged male mice. Bone 133:115255. doi: 10.1016/j.bone.2020.115255

PubMed Abstract | CrossRef Full Text | Google Scholar

Galea, G. L., Hannuna, S., Meakin, L. B., Delisser, P. J., Lanyon, L. E., and Price, J. S. (2015). Quantification of alterations in cortical bone geometry using site specificity software in mouse models of aging and the responses to ovariectomy and altered loading. Front. Endocrinol. (Lausanne) 6:52. doi: 10.3389/fendo.2015.00052

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, J. M., Ayati, B. P., Ramakrishnan, P. S., and Martin, J. A. (2012). Towards a new spatial representation of bone remodeling. Math. Biosci. Eng. 9, 281–295. doi: 10.3934/mbe.2012.9.281

PubMed Abstract | CrossRef Full Text | Google Scholar

Halloran, B. P., Ferguson, V. L., Simske, S. J., Burghardt, A., Venton, L. L., and Majumdar, S. (2002). Changes in bone structure and mass with advancing age in the male C57BL/6J mouse. J. Bone Miner. Res. 17, 1044–1050. doi: 10.1359/jbmr.2002.17.6.1044

PubMed Abstract | CrossRef Full Text | Google Scholar

Hildebrand, T., and Rüegsegger, P. (1997). A new method for the model-independent assessment of thickness in three-dimensional images. J. Microsc. 185, 67–75. doi: 10.1046/j.1365-2818.1997.1340694.x

CrossRef Full Text | Google Scholar

Hjelmstad, K. D. (2005). Fundamentals of Structural Mechanics. Boston, MA: Springer, doi: 10.1007/b101129

CrossRef Full Text | Google Scholar

Huiskes, R., Weinans, H., Grootenboer, H. J., Dalstra, M., Fudala, B., and Slooff, T. J. (1987). Adaptive bone-remodeling theory applied to prosthetic-design analysis. J. Biomech. 20, 1135–1150. doi: 10.1016/0021-9290(87)90030-3

CrossRef Full Text | Google Scholar

Javaheri, B., Bravenboer, N., Bakker, A. D., van der Veen, A., de Souza, R. L., Saxon, L., et al. (2019). “In Vivo Models of Mechanical Loading BT,” in Bone Research Protocols, ed. A. I. Idris (New York, NY: Springer New York), 369–390. doi: 10.1007/978-1-4939-8997-3_22

CrossRef Full Text | Google Scholar

Kohles, S. S., Bowers, J. R., Vailas, A. C., and Vanderby, R. (1997). Ultrasonic wave velocity measurement in small polymeric and cortical bone specimens. J. Biomech. Eng. 119, 232–236. doi: 10.1115/1.2796085

CrossRef Full Text | Google Scholar

Lambers, F. M., Kuhn, G., Weigt, C., Koch, K. M., Schulte, F. A., and Müller, R. (2015). Bone adaptation to cyclic loading in murine caudal vertebrae is maintained with age and directly correlated to the local micromechanical environment. J. Biomech. 48, 1179–1187. doi: 10.1016/j.jbiomech.2014.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Lerebours, C., Buenzli, P. R., Scheiner, S., and Pivonka, P. (2016). A multiscale mechanobiological model of bone remodelling predicts site-specific bone loss in the femur during osteoporosis and mechanical disuse. Biomech. Model. Mechanobiol. 15, 43–67. doi: 10.1007/s10237-015-0705-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y., Boudiffa, M., Dall’Ara, E., Bellantuono, I., and Viceconti, M. (2016). Development of a protocol to quantify local bone adaptation over space and time: quantification of reproducibility. J. Biomech. 49, 2095–2099. doi: 10.1016/j.jbiomech.2016.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

McClung, M. R. (2016). Cancel the denosumab holiday. Osteoporos. Int. 27, 1677–1682. doi: 10.1007/s00198-016-3553-3

PubMed Abstract | CrossRef Full Text | Google Scholar

McKenzie, J. A., and Silva, M. J. (2011). Comparing histological, vascular and molecular responses associated with woven and lamellar bone formation induced by mechanical loading in the rat ulna. Bone 48, 250–258. doi: 10.1016/j.bone.2010.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Meakin, L. B., Price, J. S., and Lanyon, L. E. (2014). The contribution of experimental in vivo models to understanding the mechanisms of adaptation to mechanical loading in bone. Front. Endocrinol. (Lausanne) 5:154. doi: 10.3389/fendo.2014.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

Minisola, S., Cipriani, C., Della Grotta, G., Colangelo, L., Occhiuto, M., Biondi, P., et al. (2019). Update on the safety and efficacy of teriparatide in the treatment of osteoporosis. Ther. Adv. Musculoskelet. Dis. 11:1759720X19877994. doi: 10.1177/1759720X19877994

PubMed Abstract | CrossRef Full Text | Google Scholar

Moustafa, A., Sugiyama, T., Prasad, J., Zaman, G., Gross, T. S., Lanyon, L. E., et al. (2012). Mechanical loading-related changes in osteocyte sclerostin expression in mice are more closely associated with the subsequent osteogenic response than the peak strains engendered. Osteoporos. Int. 23, 1225–1234. doi: 10.1007/s00198-011-1656-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ozcivici, E., Luu, Y. K., Adler, B., Qin, Y. X., Rubin, J., Judex, S., et al. (2010). Mechanical signals as anabolic agents in bone. Nat. Rev. Rheumatol. 6, 50–59. doi: 10.1038/nrrheum.2009.239

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereira, A. F., Javaheri, B., Pitsillides, A. A., and Shefelbine, S. J. (2015). Predicting cortical bone adaptation to axial loading in the mouse tibia. J. R. Soc. Interf. 12:0590. doi: 10.1098/rsif.2015.0590

PubMed Abstract | CrossRef Full Text | Google Scholar

Pickering, E., Silva, M. J., Delisser, P., Brodt, M. D., Gu, Y. T., and Pivonka, P. (2021). Estimation of load conditions and strain distribution for in vivo murine tibia compression loading using experimentally informed finite element models. J. Biomech. 115:110140. doi: 10.1016/j.jbiomech.2020.110140

PubMed Abstract | CrossRef Full Text | Google Scholar

Pivonka, P., Park, A., and Forwood, M. (2018). “Functional adaptation of bone: the mechanostat and beyond,” in CISM International Centre for Mechanical Sciences, Courses and Lectures, ed. P. Pivonka (Cham: Springer), 1–60. doi: 10.1007/978-3-319-58845-2_1

CrossRef Full Text | Google Scholar

Razi, H., Birkhold, A. I., Zaslansky, P., Weinkamer, R., Duda, G. N., Willie, B. M., et al. (2015). Skeletal maturity leads to a reduction in the strain magnitudes induced within the bone: a murine tibia study. Acta Biomater. 13, 301–310. doi: 10.1016/j.actbio.2014.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, B. C., Arredondo Carrera, H. M., Zanjani-pour, S., Boudiffa, M., Wang, N., Gartland, A., et al. (2020). PTH(1–34) treatment and/or mechanical loading have different osteogenic effects on the trabecular and cortical bone in the ovariectomized C57BL/6 mouse. Sci. Rep. 10:8889. doi: 10.1038/s41598-020-65921-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Robling, A. G., Burr, D. B., and Turner, C. H. (2001). Recovery periods restore mechanosensitivity to dynamically loaded bone. J. Exp. Biol. 204, 3389–3399.

Google Scholar

Robling, A. G., Castillo, A. B., and Turner, C. H. (2006). Biomechanical and molecular regulation of bone remodeling. Annu. Rev. Biomed. Eng. 8, 455–498. doi: 10.1146/annurev.bioeng.8.061505.095721

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubin, C. T., and Lanyon, L. E. (1984). Regulation of bone formation by applied dynamic loads. J. Bone Joint Surg. Am. 66, 397–402.

Google Scholar

Rubin, C. T., and Lanyon, L. E. (1985). Regulation of bone mass by mechanical strain magnitude. Calcif. Tissue Int. 37, 411–417.

Google Scholar

Scheuren, A. C., Vallaster, P., Kuhn, G. A., Paul, G. R., Malhotra, A., Kameo, Y., et al. (2020). Mechano-regulation of trabecular bone adaptation is controlled by the local in vivo environment and logarithmically dependent on loading frequency. Front. Bioeng. Biotechnol. 8:566346. doi: 10.3389/fbioe.2020.566346

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulte, F. A., Ruffoni, D., Lambers, F. M., Christen, D., Webster, D. J., Kuhn, G., et al. (2013). Local mechanical stimuli regulate bone formation and resorption in mice at the tissue level. PLoS One 8:e62172. doi: 10.1371/journal.pone.0062172

PubMed Abstract | CrossRef Full Text | Google Scholar

Srinivasan, S., Weimer, D. A., Agans, S. C., Bain, S. D., and Gross, T. S. (2002). Low-magnitude mechanical loading becomes osteogenic when rest is inserted between each load cycle. J. Bone Miner. Res. 17, 1613–1620. doi: 10.1359/jbmr.2002.17.9.1613

PubMed Abstract | CrossRef Full Text | Google Scholar

Stadelmann, V. A., Bonnet, N., and Pioletti, D. P. (2011). Combined effects of zoledronate and mechanical stimulation on bone adaptation in an axially loaded mouse tibia. Clin. Biomech. 26, 101–105. doi: 10.1016/j.clinbiomech.2010.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugiyama, T., Meakin, L., Browne, W., Galea, G., Price, J., and Lanyon, L. (2012). Bones’ adaptive response to mechanical loading is essentially linear between the low strains associated with disuse and the high strains associated with the lamellar/woven bone transition. J. Bone Miner. Res. 27, 1784–1793. doi: 10.1002/jbmr.1599

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugiyama, T., Meakin, L. B., Galea, G. L., Jackson, B. F., Lanyon, L. E., Ebetino, F. H., et al. (2011). Risedronate does not reduce mechanical loading-related increases in cortical and trabecular bone mass in mice. Bone 49, 133–139. doi: 10.1016/j.bone.2011.03.775

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugiyama, T., Price, J. S., and Lanyon, L. E. (2010). Functional adaptation to mechanical loading in both cortical and cancellous bone is controlled locally and is confined to the loaded bones. Bone 46, 314–321. doi: 10.1016/j.bone.2009.08.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Tiwari, A. K., Kumar, R., Tripathi, D., and Badhyal, S. (2018). In silico modeling of bone adaptation to rest-inserted loading: strain energy density versus fluid flow as stimulus. J. Theor. Biol. 446, 110–127. doi: 10.1016/j.jtbi.2018.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Trichilo, S. (2018). Assessment of the Anabolic Effects of PTH Drug Treatment and Mechanical Loading on Bone Using High-Resolution Imaging and in Silico Modelling. Ph. D. thesis. Melbourne VIC: University of Melbourne.

Google Scholar

Webster, D., Wirth, A., Van Lenthe, G. H., and Müller, R. (2012). Experimental and finite element analysis of the mouse caudal vertebrae loading model: prediction of cortical and trabecular bone adaptation. Biomech. Model. Mechanobiol. 11, 221–230. doi: 10.1007/s10237-011-0305-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolff, J. (1893). Das Gesetz der Transformation der Knochen. Berlin: Pro Business, doi: 10.1055/s-0028-1144106

CrossRef Full Text | Google Scholar

Zhang, Y., Dall’Ara, E., Viceconti, M., and Kadirkamanathan, V. (2019). A new method to monitor bone geometry changes at different spatial scales in the longitudinal in vivo μCT studies of mice bones. PLoS One 14:e0219404. doi: 10.1371/journal.pone.0219404

PubMed Abstract | CrossRef Full Text | Google Scholar