Predicting cortical bone adaptation to axial loading in the mouse tibia

Abstract

The development of predictive mathematical models can contribute to a deeper understanding of the specific stages of bone mechanobiology and the process by which bone adapts to mechanical forces. The objective of this work was to predict, with spatial accuracy, cortical bone adaptation to mechanical load, in order to better understand the mechanical cues that might be driving adaptation. The axial tibial loading model was used to trigger cortical bone adaptation in C57BL/6 mice and provide relevant biological and biomechanical information. A method for mapping cortical thickness in the mouse tibia diaphysis was developed, allowing for a thorough spatial description of where bone adaptation occurs. Poroelastic finite-element (FE) models were used to determine the structural response of the tibia upon axial loading and interstitial fluid velocity as the mechanical stimulus. FE models were coupled with mechanobiological governing equations, which accounted for non-static loads and assumed that bone responds instantly to local mechanical cues in an on-off manner. The presented formulation was able to simulate the areas of adaptation and accurately reproduce the distributions of cortical thickening observed in the experimental data with a statistically significant positive correlation (Kendall's τ rank coefficient τ = 0.51, p < 0.001). This work demonstrates that computational models can spatially predict cortical bone mechanoadaptation to a time variant stimulus. Such models could be used in the design of more efficient loading protocols and drug therapies that target the relevant physiological mechanisms.

Keywords: bone mechanobiology; cortical thickness; fluid flow; functional adaptation; mouse tibia.

Figure 1.

(a) The loading model used for this study: axial compression of the murine tibia (adapted from [23]). (b) Non-dimensional coordinate, Z, corresponds to the axial space between the proximal (Z = 0) and distal (Z = 1) tibia–fibula junctions. (c–e) Cross-sectional contours of left (non-loaded, grey line) and right (loaded, black line) tibiae for a particular specimen. Anatomical landmarks: interosseous crest (IC), proximal tibial crest (PTC), soleal line (SL), tibial ridge (TR). (c) Z = 0.2, (d) Z = 0.5 and (e) Z = 0.8.

Figure 2.

Steps of script for the calculation of cortical thickness and endosteal distance from the centroid. (a) Initial cross section, (b) peri- and endosteal boundaries and (c) map of distances from the endosteum (units in micrometres).

Figure 3.

Angular (θ) and longitudinal (Z) coordinates used in the cylindrical coordinates system.

Figure 4.

Different (mechanical stimulus versus adaptation) relationship curves considered in this study: (a) generic trilinear curve, (b) trilinear curve capped with maximum remodelling rates of apposition and resorption, (c) on–off relationship. The apposition and resorption limits of the homeostatic interval, Ψ_A and Ψ_R, dictate the tissue response to mechanical stimulus Ψ. The step curve in (c) obtained the most accurate predictions of cortical adaptation.

Figure 5.

Cortical bone adaptation model. An externally applied force elicits a mechanical signal/stimulus, Ψ, in the bone. If this stimulus is above a certain apposition threshold, Ψ_A, then an adaptive response is generated (Λ = 1). This response is integrated over time to give the adaptation rate (u̇).

Figure 6.

Adaptation algorithm flowchart. Shade indicative of software used: Mimics (green), ABAQUS (blue) and Matlab (red). Steps inside dashed line are automated.

Figure 7.

FE model used to determine the mechanical environment. (a) Medial view of the material sections considered: poroelastic cortical bone (off-white), elastic cortical bone (blue) and fibular growth plate (red). (b) Lateral view of the boundary conditions applied and vertical alignment of the model assembly with the z-direction. (c) Medial view of the region where cortical bone adaptation was considered, demarcated by the cloud of nodes highlighted in red.

Figure 8.

Contour plots of the longitudinal strains calculated in the FE models and DIC measurements in eight-week-old mice of the same strain and gender by Sztefek et al. [34]. FEA used the geometry of a non-adapted tibia under an axial peak load of 12 N. (DIC images adapted from [34] with permission from Elsevier.)

Figure 9.

Contour plots of fluid velocity, V, obtained at t = 0.025 s (peak load) for different Z positions in the initial geometry. Higher fluid velocity magnitudes obtained in the medial surface and around the interosseous crest. (a) Z = 0.1, (b) Z = 0.3, (c) Z = 0.5, (d) Z = 0.7 and (e) Z = 0.9.

Figure 10.

(a–c) Cross-sectional contours of left (non-loaded, grey line) and right (loaded, black line) tibiae, and respective polar representation cortical thickness, Th, for (d–f) the in vivo data measurements and (g–i) predictions using Ψ_A = 50 MPa. The red asterisk (c and i) is placed at Z = 0.7 and θ = 135° and highlights a region where in silico adaptation was overestimated (Th units in micrometres). (a) Z = 0.3, (b) Z = 0.5, (c) Z = 0.7; (d) micro-CT scan, Z = 0.3, (e) micro-CT scan, Z = 0.5, (f) micro-CT scan, Z = 0.7; (g) simulation, Z = 0.3, (h) simulation, Z = 0.5 and (i) simulation, Z = 0.7.

Figure 11.

Colour map of changes in cortical thickness, ΔTh/Th, measured for in vivo data (a) against the respective predictions using Ψ_A = 50 MPa (b). The red asterisk is placed at Z = 0.7 and θ = 135°. Artefacts in the experimental data calculations predominate in the region masked over, outside the yellow dashed line. (a) micro-CT scan and (b) simulation.

Figure 12.

Kendall's τ rank correlation coefficient calculated between in vivo and in silico ΔTh/Th measurements obtained for the values of Ψ_A considered in this study. Simulations with an apposition threshold value of Ψ_A = 3 × 10⁻² ms⁻¹ obtained the strongest positive correlation (p < 0.001).