Mineral density differences between femoral cortical bone and trabecular bone are not explained by turnover rate alone.

Bone mineral density distributions (BMDDs) are a measurable property of bone tissues that depends strongly on bone remodelling and mineralisation processes. These processes can vary significantly in health and disease and across skeletal sites, so there is high interest in analysing these processes from experimental BMDDs. Here, we propose a rigorous hypothesis-testing approach based on a mathematical model of mineral heterogeneity in bone due to remodelling and mineralisation, to help explain differences observed between the BMDD of human femoral cortical bone and the BMDD of human trabecular bone. Recent BMDD measurements show that femoral cortical bone possesses a higher bone mineral density, but a similar mineral heterogeneity around the mean compared to trabecular bone. By combining this data with the mathematical model, we are able to test whether this difference in BMDD can be explained by (i) differences in turnover rate; (ii) differences in osteoclast resorption behaviour; and (iii) differences in mineralisation kinetics between the two bone types. We find that accounting only for differences in turnover rate is inconsistent with the fact that both BMDDs have a similar spread around the mean, and that accounting for differences in osteoclast resorption behaviour leads to biologically inconsistent bone remodelling patterns. We conclude that the kinetics of mineral accumulation in bone matrix must therefore be different in femoral cortical bone and trabecular bone. Although both cortical and trabecular bone are made up of lamellar bone, the different mineralisation kinetics in the two types of bone point towards more profound structural differences than usually assumed.

Introduction

The mineral heterogeneity of bone tissues is the result of remodelling and mineralisation processes, which are known to vary significantly in health and disease and across skeletal sites (Boivin and Meunier, 2002b; Fratzl et al., 2007; Kingsmill et al., 2007; Roschger et al., 2008; Misof et al., 2014b). Bone remodelling periodically replaces old bone with an unmineralised collagen-rich matrix, which is subsequently infiltrated with a mineral phase. The growth of mineral crystals in this phase gradually confers to bone its stiffness and strength (Parfitt, 1983; Currey, 1999). Regions of bone formed at different times achieve various degrees of mineralisation and are visualised as different grey-level intensities in bone scans obtained by X-ray absorption or quantitative backscattered electron imaging (qBEI) (Roschger et al., 1998; Boivin and Meunier, 2002b; Busse et al., 2009; Nuzzo et al., 2002; Rohrbach et al., 2012; Sansalone et al., 2010; Dierolf et al., 2010). The heterogeneity of mineral density in a bone scan gives an indication of bone's renewal history and provides an indirect measure of remodelling and mineralisation processes (Buenzli, 2016). Bone mineral density distributions (BMDDs) quantify this mineral heterogeneity experimentally as frequency distributions of calcium content constructed from qBEI scans (Roschger et al., 2003, Roschger et al., 2008).

In adult trabecular bone, the BMDD is independent of age, ethnicity, sex, and skeletal site (transiliac bone, vertebrae, femoral neck, femoral head, and patella) (Roschger et al., 2003, Roschger et al., 2008), and the relatively low inter-individual variation has enabled the definition of a reference trabecular BMDD of healthy adults (Roschger et al., 2003). In several diseases, the BMDD deviates from this reference because of differences in remodelling and/or mineralisation processes (Roschger et al.,2008). In cortical bone, the mean mineral content depends on skeletal site (Kingsmill et al., 2007), so that cortical BMDDs may differ from the reference trabecular BMDD in health too.

Novel measurements of cortical BMDDs from human femur midshafts in healthy adults reveal a higher degree of mineralisation than the reference trabecular BMDD, but a similar spread of mineral heterogeneity around the most frequently occurring calcium content (position of the peak of the BMDD) (Buenzli et al., 2018). To understand the reason for the difference in mineral distribution between trabecular bone and femoral cortical bone in healthy adults, we propose in this paper a systematic hypothesis-testing analysis of BMDDs using an established mathematical model of bone mineral heterogeneity due to remodelling and mineralisation. This mathematical model allows us to quantify BMDD signatures in terms of descriptive parameters of remodelling and mineralisation processes (Ruffoni et al., 2007, Ruffoni et al., 2008).

Shifts of BMDDs and of average bone mineral densities towards higher and lower mineral densities have previously been attributed to lower and higher birth rates of basic multicellular units (BMUs), respectively (Boivin and Meunier, 2002a; Misof et al., 2012), and acceleration and deceleration of mineral accumulation in specific bone diseases (Nawrot-Wawrzyniak et al., 2009). Mineralisation kinetics, i.e. the accumulation of mineral density with time in a small volume of newly formed bone, has been reported to be similar in cortical and trabecular bones in healthy dogs (Marotti et al., 1972), and in the iliac crest of healthy ewes (Bala et al., 2010). Since turnover rates are normally lower in femoral cortical bone than in trabecular bone (Manolagas, 2000; Parfitt, 2002), this difference could explain the shift in BMDD towards higher mineral densities seen in human femoral cortical bone. A lower turnover rate provides more time for secondary mineralisation to take place, which increases the average bone mineral density (Boivin and Meunier, 2002a). However, lower turnover rates also lead to broader BMDD peaks (Roschger et al., 2008; Ruffoni et al., 2008; Buenzli et al., 2018). For the femoral cortical BMDD to possess a spread around the mean similar to that of the trabecular BMDD, additional processes than turnover rate may need to be accounted for. We thus explore systematically and quantitatively further hypotheses related to differences in resorption patterns of remodelling, and to the timescale of mineralisation kinetics. One of the novelties of the mathematical analysis presented in this paper is to investigate whether resorption behaviour of osteoclasts targeting specific calcium content could explain the differences between the BMDDs. Resorption patterns are likely to be different in cortical bone, where resorption may reach any mineral density by tunnelling of new Haversian canals, compared to trabecular bone, where resorption predominantly reaches weakly mineralised bone at the trabeculae's surfaces (Martin et al., 1998; Eriksen, 2010).

Materials and methods

Backscattered electron imaging and experimental BMDD

Midshaft femur cross-sections were extracted post mortem at the Department of Forensic Medicine of the Medical University of Vienna from four adult women (48, 50, 55 and 56 years old), who died unexpectedly with no known bone pathologies. The study was performed in accordance with the ethic commission board of the institution (EK#: 1757/2013). A single femur cross-section was extracted from each donor. The samples were frozen at −20 °C immediately after extraction. Prior to measurements, the samples were unfrozen, dehydrated and embedded in polymethylmethacrylate (PMMA) as described elsewhere (Roschger et al., 1998). Using the scanning electron microscope (Zeiss supra 40, Oberkochen, Germany) in backscattered mode, qBEI scans of in average of 40 mm2 bone area per sample covering medial or lateral cortical bone were acquired with a pixel resolution of 1.73 μm. The imaging experimental setup was calibrated and checked for measurement reliability to convert grey intensity values into calcium weight percent (wt% Ca), and further to extract BMDDs following standard procedures (Roschger et al., 1998; Lukas et al., 2011). The evaluation was done twice, first by considering the whole cortex and then, by excluding the endosteal and periosteal regions. Obtained differences were minor and, therefore, only results from a whole cortex evaluation are presented. A mean femoral cortical BMDD (Fm.BMDD) was constructed from the very similar individual BMDDs by calculating the mean values for each histogram bin (Buenzli et al., 2018).

The reference trabecular BMDD (Tb.BMDD) used in this paper corresponds to the average of 52 BMDDs coming from different individuals, and different skeletal sites from a previous study (Roschger et al., 2003). Although the number of individuals investigated here to define (Fm.BMDD) is lower than the cohort used to establish the reference trabecular BMDD, the cortical BMDD statistics is substantially improved by the large areas measured. In Roschger et al. (2003), 5–6 regions of tissue area 5 mm2 were used per cancellous bone sample, but the analysed bone area is reduced due to the much smaller BV/TV of trabecular bone compared to femoral cortical bone. In the femur, the degree of mineralisation is known to depend slightly on the region of interest along the periosteal-endosteal axis (Sansalone et al., 2010). However, the variation in BMDD within the cross-sections was small enough to justify the consideration of whole cross-section BMDDs.

The Tb.BMDD and Fm.BMDD were smoothed by the ‘loess’ Matlab procedure (MATLAB, 2017) with a span of 1.05 wt% Ca. This smoothing enabled better numerical evaluations in the mathematical analyses by alleviating the poor signal-to-noise ratio of the BMDD at very low and very high calcium content (Lukas et al., 2011).

The trabecular and femoral cortical BMDDs exhibit an overall bell-shaped curve (Fig. 1). To quantify this shape, we calculated the following shape parameters introduced in Refs (Roschger et al., 2003, Roschger et al., 2008): CaPEAK, the calcium content at the peak of the BMDD; CaMEAN, the mean of the BMDD, which corresponds to the average calcium content of the bone sample; and CaWIDTH, the full width at half maximum of the BMDD, which characterises the heterogeneity in mineral density of the sample around the distribution peak (Roschger et al., 2003, Roschger et al., 2008). We also calculated the skewness CaSKEW and the kurtosis CaKURT of the distribution, which characterise the asymmetry and the peakedness of the distribution, respectively. Skewness and kurtosis are calculated on a restricted interval where the BMDD is larger than 1% of the peak height.

Fig. 1

Experimental BMDDs: Reference trabecular BMDD (Tb.BMDD, thick black line) with the ±SD interval (grey lines); adult femoral cortical mean BMDD (Fm.BMDD: thick orange line, with corresponding individual BMDDs as thin orange lines). The BMDDs are normalised to the evaluated bone area (100%).

The strength of the BMDD analysis presented in this study, is to address rather late stages of mineralisation, which are difficult to investigate with labelling techniques.

Mathematical model of the BMDD

To evaluate the influence of turnover rate and mineralisation kinetics on BMDDs, we use the mathematical model developed by Ruffoni et al., 2007, Ruffoni et al., 2008. This mathematical model is based on the general, mechanistic principles of the balance of elementary mineralised tissue volumes in a region of interest of bone tissue (TV) during the evolution of the tissue under remodelling and mineralisation processes. The model applies to a region of bone tissue undergoing remodelling, in which mineral density evolves by the processes of (1) new bone formation; (2) bone mineralisation; and (3) bone resorption. Our main assumption is that the mineralisation kinetics is uniform within TV, i.e., independent of space. The influence of these fundamental processes on the BMDD is summarised in Fig. 2. We present the main features of this mathematical model below and refer the reader to Appendix A for further details. In the mathematical equations, the BMDD is denoted by ρ(c), where c is the calcium content. The quantity ρ(c)dc represents the volume of bone within the sample that possesses a calcium content in the histogram bin [c, c + dc). The remodelling and mineralisation processes are fully characterised in the model by three elements: the mineralisation rate, the resorption rate, and the formation rate, as described below.

Fig. 2

(a) qBEI image of iliac crest trabecular bone. Distinct bone packets within trabeculae are seen as regions of distinct grey colours. (b) Schematic representation of the zoomed qBEI image, where the different degree of mineralisation between the bone packets are highlighted. The influence of mineralisation without remodelling, and of remodelling without mineralisation on the calcium distribution in this region of bone is schematised in (c) and (d), respectively. The corresponding changes in BMDD due to these processes are shown in (e)–(f): (e) mineralisation rate q(c) (orange arrows) pushes the BMDD towards higher calcium content; the peak width narrows due to the slowdown of mineralisation kinetics as mineralisation proceeds (Buenzli et al., 2018); and (f) bone resorption rate r(c) (green arrows) removes bone and thus lowers the BMDD while bone formation (black bar) introduces unmineralised bone at c=0.

The mineralisation rate, denoted by q(c) (in wt% Ca per unit time), describes how fast calcium accumulates in a microscopic region of bone of calcium content c. This function is defined such that during a small time increment Δt, calcium content increases by Δc = q(c)Δt. The mineralisation rate q(c) entirely determines the kinetics of mineralisation c(t), i.e. how the calcium content of this region accumulates with time. Conversely, the mineralisation kinetics c(t) entirely determines the mineralisation rate q(c), see Eq. (A.8) (Ruffoni et al., 2007). The effect of mineralisation, which occurs at different speeds q(c) within a sample because of the heterogeneity of mineral density, is to push the BMDD towards higher calcium content, as illustrated in Fig. 2(e).

The resorption rate r(c) is the probability per unit time that bone matrix with calcium content c is resorbed by osteoclasts. Different osteoclast resorption patterns can be captured by modifying this function. The effect of this function in the mathematical model is to lower the height of the BMDD with different propensities depending on c, as illustrated in Fig. 2(f). A calcium-independent resorption rate represents random resorption, i.e. each region of bone is as likely to be resorbed, irrespective of its calcium content. The quantity r(c)ρ(c)dc is the volume of bone with calcium content in [c, c + dc) resorbed per unit time. Summing r(c)ρ(c)dc over all the histogram bins gives the total volume of bone resorbed per unit time in the sample, see Eq. (A.4).

Finally, the formation rate is the total volume of bone formed per unit time by the osteoblasts. Since osteoblasts form unmineralised matrix, the calcium content deposited by the osteoblasts is taken to be 0 wt% Ca.

Because we investigate differences in BMDDs in healthy adults, we assume that bone resorption and bone formation are balanced. The bone volume formed per unit time and the bone volume resorbed per unit time are taken to be identical, and we express these quantities by the turnover rate χ, defined as the fraction of bone volume (BV) renewed per unit time (Eq. (A.4)) (Eriksen, 1986). In humans, turnover rate depends on bone type and skeletal site. In this study, we take the turnover rate in trabecular bone to be χTb=20 % /year, and the turnover rate in femoral cortical bone to be χFm=5 % /year (Eriksen et al., 2002; Parfitt, 1983, Parfitt, 2002). The functions and parameters used in this model are listed in Table 1.

Table 1

Nomenclature.

Symbol	Description	Value/Unit
c	Calcium content	wt% Ca
ρ(c)	Bone Mineral Density Distribution (Tb.BMDD; Fm.BMDD)	mm3/wt% Ca
q(c)	Mineralisation rate	wt% Ca/year
r(c)	Resorption function	year−1
c(t)	Mineralisation kinetics	wt% Ca
χTb	Turnover rate in trabecular bone	20%/year
χFm	Turnover rate in femoral cortical bone	5%/year

The conservation of bone volume leads to a steady state of the BMDD, due to an equilibrium between (i) the continuous generation of new unmineralised bone due to formation; (ii) the continuous drift of bone towards higher calcium content due to mineralisation; and (iii) the continuous lowering of the BMDD due to resorption (Fig. 2). In steady state, the mathematical model expresses a relationship between the BMDD ρ(c), the mineralisation rate q(c), the resorption rate r(c) and the turnover rate χ (Ruffoni et al., 2007). Both mineralisation rate and resorption rate combine to make up a BMDD signature (see Eq. (A.5)). To extract the mineralisation rate q(c) from measurements of the BMDD (Eq. (A.6)), assumptions need to be made on the resorption rate r(c). Equivalently, to extract the resorption rate r(c) from measurements of the BMDD (Eq. (A.7)), assumptions need to be made on the mineralisation rate q(c). To explain the differences observed between Fm.BMDD and Tb.BMDD, we thus test a series of assumptions on differences that may exist in turnover rate, resorption rate, and mineralisation rate between these two bone types.

Results

The qBEI measurements show that the mean femoral cortical BMDD (Fm.BMDD) is of similar shape as the reference trabecular BMDD of healthy adults (Tb.BMDD), but it is shifted towards higher calcium contents (Fig. 1). This observation is described quantitatively by our measurements of peak shape and location. Table 2 shows that the full width at half maximum CaWIDTH, skewness CaSKEW, and kurtosis CaKURT are nearly identical for Fm.BMDD and Tb.BMDD. The difference in value of these parameters is well within the standard deviations of the measurements. However, the peak location CaPEAK of Fm.BMDD is shifted by about 3 wt% Ca, well beyond the standard deviations of CaPEAK.

Table 2

Shape parameters of the femoral cortical BMDD from adults, the reference trabecular BMDD, and their difference. The trabecular values are taken from Ref. (Roschger et al., 2008).

BMDD -parameter	Cortical	Trabecular	Difference
	Mean values	Mean values	(Cortical − Trabecular)
	±StD	±StD
CaPEAK [wt% Ca]	25.95 ± 0.22	22.94 ± 0.39	+3.01
CaMEAN [wt% Ca]	25.05 ± 0.28	22.20 ± 0.45	+2.85
CaWIDTH [wt% Ca]	3.20 ± 0.10	3.35 ± 0.34	−0.15
CaSKEW [−]	−0.91± 0.21	−0.90± 0.41	−0.01
CaKURT [−]	5.16 ± 0.68	5.17 ± 1.31	−0.01

In the following, we explore computationally a series of hypotheses on the intrinsic remodelling and mineralisation processes of femoral cortical and trabecular bone with the aim to understand how the mineral heterogeneity in these bone types can be distributed similarly, but around different peak positions.

The difference between femoral cortical and reference trabecular BMDDs is only due to the difference in turnover rate between these two bone types.

The motivation for this hypothesis is the well-known difference in turnover between bone types (Parfitt, 2002), and the suspected influence that turnover exerts on the BMDD (Boivin and Meunier, 2002a; Ruffoni et al., 2008). As a first step, we therefore assume that.

Resorption rate is slower in cortical bone than in trabecular bone due to the difference in turnover rate, i.e., the total volume fraction resorbed per unit time is lower in cortical bone;

Resorption rate is unbiased in the two bone types, i.e., osteoclasts resorb bone independently of calcium content (random resorption);

Mineralisation rate is identical in the two bone types.

To test Hypothesis 1, we check whether Fm.BMDD and Tb.BMDD measured experimentally are compatible with BMDDs calculated from Eq. (A.5) with the three assumptions above. When resorption is random, the resorption rate functions r(c) are constant and equal to the turnover rates χFm or χTb, respectively (see Eqs. (A.3), (A.4)). Two possible choices can then be made to determine the common mineralisation rate q(c) from the experimental data. The first choice corresponds to the mineralisation rate determined from Tb.BMDD and χTb using Eq. (A.6), as was done in Ref. (Ruffoni et al., 2007). The second choice corresponds to the mineralisation rate determined from Fm.BMDD and χFm using Eq. (A.6). By means of Eq. (A.5), Hypothesis 1 then allows us to estimate computationally Fm.BMDD from the experimental Tb.BMDD, denoted by “Fm.BMDD (Hyp.1)” in Fig. 3, and to estimate computationally Tb.BMDD from the experimental Fm.BMDD, denoted by “Tb.BMDD (Hyp.1)” in Fig. 3. Fig. 3 shows that the BMDDs computed using either of these choices of mineralisation rates differ significantly from the BMDDs measured experimentally. In particular, Fm.BMDD (Hyp.1) has a CaPEAK of 24.5 wt% Ca and a CaWIDTH of 5.15 wt% Ca, and Tb.BMDD (Hyp.1) has a CaPEAK of 24.5 wt% Ca and a CaWIDTH of 2.80 wt% Ca (compare with the experimental values in Table 2).

Fig. 3

Comparison between experimental BMDDs (thick lines) and BMDDs computed with Hypothesis 1 (thin lines). There is a large discrepancy between the measured femoral cortical BMDD (thick orange) and the BMDD computed with cortical turnover rate and mineralisation rate extracted from the reference trabecular BMDD (thin orange). There is also a large discrepancy between the reference trabecular BMDD (thick black) and the BMDD computed with trabecular turnover rate and mineralisation rate extracted from the femoral cortical BMDD (thin black).

The mathematical model, therefore, shows that turnover rate influences CaPEAK as expected, i.e., CaPEAK is low at high turnover rate and CaPEAK is high at low turnover rate. But turnover rate has also a strong influence on peak width CaWIDTH. Both experimental BMDDs have a similar CaWIDTH (Table 2), so that no value of turnover rate difference between femoral cortical bone and trabecular bone is able to match jointly CaPEAK and CaWIDTH of the computed BMDDs with those of the measured BMDDs (Ruffoni et al., 2007; Buenzli et al., 2018). Thus, the difference in mineral distribution between the two bone types cannot be attributed to a change in turnover rate alone.

The difference between femoral cortical and trabecular BMDDs is due to a difference in turnover rate, and to an osteoclastic bone resorption biased towards specific calcium contents in one of the two bone types.

The motivation for this hypothesis is the likely difference in resorption pattern between femoral cortical bone and trabecular bone. As a second step in our computational hypothesis testing, we still assume random resorption in trabecular bone, but explore whether there are biased resorption patterns in femoral cortical bone such that the BMDD computed from Eq. (A.5) matches the measured femoral cortical BMDD. Biased resorption towards specific calcium contents may represent a combination of mineral-density-dependent speed of resorption, and preferential initiation of resorption in regions or surfaces of bone with specific mineral densities. Mathematically, these situations are captured by letting the resorption rate be a general function r(c) of calcium content. The unique resorption function r(c) that is consistent with the measured femoral cortical BMDD can be computed using Eq. (A.7). Fig. 3 shows that this function is negative for a wide range of values of calcium content. Conversely, assuming random resorption in cortical bone, the unique trabecular bone resorption function such that the BMDD computed from Eq. (A.5) matches the measured reference trabecular BMDD is also negative for a range of values of calcium content (Fig. 4). Negative parts of r(c) for c>0 represent the formation of readily mineralised bone. The area under the negative region of the graphs of r(c)ρ(c)/χ shown in Fig. 4 represents the total volume of bone that would be formed readily mineralised per turnover time. Under Hypothesis 2, a significant amount of bone would be formed readily mineralised, which clearly contradicts what is known about osteoblast action and bone mineralisation. Hypothesis 2 is thereby falsified.

Fig. 4

Distribution of volume of bone resorbed per unit time under Hypothesis 2 with random resorption assumed in trabecular bone (orange), and random resorption assumed in femoral cortical bone (black). The unique resorption rate functions r(c) for which Hypothesis 2 holds are negative in a region of calcium content, which means Hypothesis 2 would only hold if there is formation by osteoblasts of readily mineralised bone.

The difference between femoral cortical and trabecular BMDDs is due to a difference in turnover rate, and to osteoclastic bone resorption biased towards specific calcium content differently in both bone types.

As a third step, we now relax the assumption that random resorption occurs in one of the bone types and we attempt to determine whether there are two general resorption rate functions rFm(c) and rTb(c), such that the BMDDs computed from Eq. (A.5) match the experimental BMDDs. Mineralisation kinetics is still assumed to be identical in both bone types as in Hypotheses 1 and 2, i.e., we still consider Assumption (3) to hold. The analysis is now complicated by the fact that mineralisation rates are determined from the experimental BMDDs under a choice of resorption functions by Eq. (A.6). The procedure employed to test Hypothesis 3, is therefore (i) to first choose a resorption function rTb(c); (ii) to deduce the mineralisation rate from Eq. (A.6) using the femoral cortical BMDD and the choice of resorption function rTb(c); and (iii) to then determine rFm(c) from Eq. (A.7) using the reference trabecular BMDD. We can also proceed conversely by choosing first rFm(c), and calculate how this choice determines rTb(c). This procedure is repeated to explore the set of all possible choices of positive resorption functions rTb(c) that lead to positive resorption functions rFm(c), and vice versa. Such a pair of positive resorption functions rTb(c), rFm(c) represents targeted osteoclast resorption behaviours that lead to the experimental BMDDs while having identical mineralisation kinetics in femoral cortical bone and trabecular bone. Below, we show that this set is strongly restricted by the requirement that both resorption functions must be positive. We find that only cortical resorption functions that strongly target very lowly mineralised bone satisfy this requirement.

The detailed mathematical analysis is more involved than in the previous cases, since it requires the exploration of a large space of possible resorption functions. The reader is referred to Appendix B for the mathematical details. To help us enforce the positivity of rFm(c) given a choice of rTb(c), we introduce an auxiliary function GTb(c) defined in Eq. (B.4). Choosing an initial resorption function rTb(c) is equivalent to choosing an auxiliary function GTb(c). The advantage is that to explore the set of all resorption functions rTb(c) that give rise to positive resorption functions rFm(c), it is sufficient to explore the set of all auxiliary functions GTb(c) below an explicit upper limit determined by the experimental BMDDs (Eq. (B.7)). Conversely, to explore the set of all resorption functions rFm(c) that give rise to positive resorption functions rTb(c), it is sufficient to explore the set of all auxiliary functions GFm(c) below an explicit upper bound.

To illustrate the difficulty of enforcing the positivity of the resorption functions, Fig. 5a, c show a number of different initial choices of rFm(c) that correspond to osteoclasts resorbing cortical bone that is preferentially (i) lowly mineralised, by assuming that r(c) = r0 + λ(e − 1) with negative exponents γ = −10, −30 and r0=5.72×106, 15, respectively; (ii) highly mineralised by assuming that r(c) = r0 + λ(e − 1) with positive exponents γ=5, 20 and r0=0.025; and (iii) of calcium content restricted to Gaussians of different centres M and widths σ, i.e., . Since the resorption functions must satisfy the turnover condition (A.4), they only depend on two independent parameters. Fig. 5b, d show that choosing functions rFm(c) that correspond to preferential resorption targeting intermediate or highly mineralised cortical bone (γ=5, 20, and M=25, 28) do not lead to trabecular resorption rate functions rTb(c) that are positive. The auxiliary functions GFm(c) corresponding to these choices of rFm(c) lie in part above the upper limit (Fig. 5b, d). The same holds for a mild resorption preference towards low calcium content (γ = −10). The only resorption functions rFm(c) that lead to positive trabecular resorption functions rTb(c) are the Gaussian centred around lowly mineralised bone (M=15), and the exponential resorbing only lowly mineralised bone (γ = −30).

Fig. 5

Resorption functions used to test Hypothesis 3; (a) exponential resorption functions; and (c) Gaussian resorption functions. The inset in (c) shows a broader view of the Gaussian centred at 15 wt% Ca. The auxiliary functions GFm(c) corresponding to these exponential and Gaussian resorption functions are shown in (b) and (d), respectively. Auxiliary functions GFm(c) are subject to the condition (B.7) derived in Appendix B, i.e., they must sit below the shaded area bounded by the dashed black line to ensure that resorption functions are positive, i.e., that no bone is formed readily mineralised. Only resorption functions strongly targeting lowly mineralised bone satisfy this constraint.

The procedure to check Hypothesis 3 thoroughly is to explore the full set of initial choices GTb(c) below the upper limit and inspect the cortical resorption functions rFm(c) determined by these choices (see Appendix B). For illustrative purposes, we show in Fig. 6 the distribution of cortical bone volume resorbed per unit time rFm(c)ρFm(c) when GTb(c) is set as the upper limit (borderline case, Eq. (B.7)). In this case too, cortical bone would be preferentially resorbed at extremely low calcium content (< 22 wt% Ca).

Fig. 6

Distribution of volume of bone resorbed per unit time in cortical bone, rFm(c)ρFm(c), obtained for certain choices of trabecular auxiliary functions GTb(c), including for when GTb(c) is set as the upper bound of the condition (B.7). The measured femoral cortical BMDD (rescaled by χFm) is also shown to emphasise that the cortical resorption behaviours obtained in these cases renew a very small volume of lowly mineralised cortical bone at an extremely high rate.

In Appendix B.1, we show rigorously that under Hypothesis 3, for any choice of auxiliary function, at least 25% of the femoral cortical bone resorbed per unit time would have a calcium content less than 21 wt% Ca, and at least 15% would have a calcium content less than 20 wt% Ca (see Eq. (B.11) and Fig. B.10). In other words, more than a quarter of all the cortical bone renewed would be targeting new, poorly mineralised bone, to the extent that more than 15% of all cortical bone would be resorbed before undergoing secondary mineralisation. These predictions clearly conflict with the concept that bone remodelling is a repair mechanism that enables bones to maintain structural integrity despite fatigue loading (Martin, 2002; Burr, 2002; Bentolila et al., 1998), especially in mechanically stimulated bone such as the femoral cortex. For all the above reasons, we therefore deem Hypothesis 3 to be falsified.

Fig. B.10

Under Hypothesis 3, the fraction α(c∗) of total volume of bone resorbed per unit time that has a calcium content less than c∗ must be greater than the lower bound in Eq. (B.11) plotted in this figure. This shows that α(21 wt % Ca) ≥ 25% and α(20 wt % Ca) ≥ 20%.

The difference between femoral cortical and trabecular BMDDs is due to a difference in turnover rate, and a difference in mineralisation kinetics between the two bone types.

Having exhausted all the possibilities in which mineralisation kinetics was assumed identical between trabecular and femoral cortical bone, we now hypothesise that mineralisation kinetics may be different in trabecular bone and femoral cortical bone, i.e., we no longer consider our assumption (3) to hold. To explore this hypothesis computationally, we revert to the assumption of random resorption in both bone types for simplicity. In this case, the mineralisation rates q(c) and mineralisation kinetic laws c(t) can be determined in each bone types independently by Eqs. (A.6), (A.8) from the experimental BMDDs, Tb.BMDD and Fm.BMDD. The mineralisation rates and kinetic laws determined in this way are shown in Fig. 7. Both mineralisation kinetics display a fast initial increase in calcium content, followed by a slower increase in calcium content with time. However, the overall rate of mineralisation in femoral cortical bone is slower than the overall rate of mineralisation in trabecular bone. The timescale of mineralisation kinetics in either bone type depends on the assumed turnover rate in Eq. (A.6). To factor out the influence of the choice of turnover rate and compare qualitatively the mineralisation kinetics in femoral cortical bone and trabecular bone, we perform a scaling in time of the mineralisation kinetics of cortical bone by the ratio of turnover rates . Remarkably, Fig. 7 shows that the time-scaled cortical mineralisation rate and time-scaled cortical mineralisation kinetic laws have similar shape as the trabecular mineralisation rate and mineralisation kinetic law, respectively, but are shifted along the calcium content axis. Accounting for this shift additionally by including an offset κ results in a nearly perfect superposition of the transformed cortical mineralisation functions and the trabecular mineralisation functions, i.e.,

Fig. 7

(a) Mineralisation rates qTb(c), qFm(c); and (b) Mineralisation kinetic laws cTb(t), cFm(c) obtained from the experimental trabecular BMDD and femoral cortical BMDD under Hypothesis 4. The mineralisation rates q(c) and kinetic laws c(t) obtained by rescaling the timescale of femoral cortical mineralisation by the ratio of turnover rate , and by shifting the calcium axis by κ ≈ 2.83 wt% Ca are also shown. The time-rescaled and calcium-shifted femoral cortical mineralisation functions q(c) and c(t) fit the trabecular mineralisation functions qTb(c) and cTb(t) remarkably well both at low calcium densities (≲20 wt% Ca, as seen more obviously in (a)), and at high calcium densities (≳20 wt% Ca, as seen more obviously in (b)).

The goodness of match between the mineralisation kinetics of trabecular and femoral cortical bone after a scaling in time and a shift in calcium content in Eqs. (1), (2) is remarkable because it holds in the full range of experimentally accurate measurements of calcium densities (Fig. 7). It is important to note that this match does not rely on the choice of ratio of turnover rates. The same match holds for other assumed values of turnover rates, so long as these values are used consistently in the other quantities that depend on them. It can be shown mathematically from Eq. (A.5) that the reason for this goodness of match is the similarities of the shapes of the two BMDDs (Fig. 1 and Table 2). By translating the femoral cortical BMDD towards lower calcium contents, such that both peak locations coincide, the two BMDDs superimpose almost perfectly. The translated femoral cortical BMDD is within the 95% confidence interval of the trabecular measurements. A sensitivity analysis of the link between BMDD shape and mineralisation kinetics is provided in Appendix C.

Discussion

Turnover rate is well-known to vary across skeletal sites and bone types. In particular, regions of femoral cortical bone of low porosity may be turned over roughly 4 times slower than trabecular bone (Eriksen et al., 2002; Parfitt, 1983, Parfitt, 2002). This has been proposed to explain the high degree of mineralisation of cortical bone compared to trabecular bone (Boivin and Meunier, 2002a). However, our analysis clearly excludes that the difference in mineral density distribution between femoral cortical bone and trabecular bone captured in the BMDDs is an effect of turnover only (Hypothesis 1, Fig. 3). By testing rigorously general hypotheses on resorption rates, we were able to exclude other possible ways by which one may explain the difference between the two BMDDs. In particular, the concept that different resorption behaviours of osteoclasts may take place in cortical versus trabecular bone, was found to be an insufficient explanation (Hypotheses 2 and 3, Fig. 4, Fig. 5, Fig. 6). The assumption that osteoclasts may resorb preferentially bone with specific calcium content can be understood as (i) a cellular preference to target microcracks, generally situated in old, highly mineralised bone (Burr, 2002; Martin, 2002); (ii) geometrical constraints related to the fact that resorption occurs first at existing surfaces in trabecular bone. Bone tissue near the bone surfaces is less mineralised than interstitial bone on average, so resorption targeting less mineralised bone represents resorption of bone layers that remain close to the bone surface. In contrast, resorption targeting more mineralised bone represents the resorption of deeper layers of bone, such as interstitial bone, which is located further away from bone surfaces. Assigning various propensities for bone tissue resorption depending on mineral content effectively captures such fine-scale spatial dependences of bone resorption. Because we make no a priori restriction on the functional relationship between calcium content and resorption rate, we do not need to specify further how mineral density correlates with spatial location within bone to account for these effects; and (iii) mechanical regulation of remodelling, which may recruit osteoclasts at different rates depending on the mechanical properties of the bone site (Burger et al., 2003) such as the local stiffness of bone matrix, which depends in particular on mineral content. Our completely general choice of possible resorption functions r(c) thereby takes into consideration in an implicit manner a broad range of influencing factors on osteoclastic behaviour.

The fact that the mineralisation kinetic laws in Fig. 7 only match after rescaling in time suggests that mineralisation kinetics is different in trabecular and femoral cortical bone. Eqs. (1), (2) indicate that the time scale of mineralisation is coupled with the time scale of bone turnover independently of bone type (see Appendix D). Our conclusion that mineralisation may proceed at a different rate in trabecular and femoral cortical bone focuses on mid/late stages of mineralisation. In a previous study, Marotti et al. studied the mineralisation kinetics by fluorescence labelling over 47 days in the dog humerus, radius, ulna, femur, tibia diaphysis, and radius metaphysis, and found no statistical differences (Marotti et al., 1972). These findings hold for early stages of mineralisation and are thus not contradictory with our conclusions. At this early stage of mineralisation, our computational model predicts a mineralisation kinetics which is very similar for the two bone types (Fig. 7). In a more recent study, Bala et al. investigated the mineralisation kinetics in iliac crest of ewes over 30 months (Bala et al., 2010), also finding no statistical differences of mineralisation kinetics between trabecular and cortical bone. This observation emphasises the importance of the skeletal site for bone mineral heterogeneity. Our conclusions are drawn from BMDD data from femoral cortical bone. The development of the cortex in the iliac crest follows an intricate development during growth (Fratzl-Zelman et al., 2009b) that, unlike femoral bone, is likely to retain traces of its growth history in terms of large-scale spatial dependences of mineral heterogeneity. Furthermore, its cortical and trabecular compartments may be regulated similarly (Misof et al., 2014a). Our analysis relies on a steady state assumption where traces of growth history have been erased by remodelling. The mathematical model also relies on the continuum modelling assumption that the tissue volume TV is large enough so that the BMDD averages many mineral heterogeneities, but small enough that there are no large-scale spatial dependences of the BMDD within TV.

The experimental basis of our study is that the peak position of Tb.BMDD and Fm.BMDD is different, but the mineral heterogeneity is similar, as quantified by CaWIDTH (Table 2). One major intent of the current work is to emphasise that this similar heterogeneity is unexpected from our current understanding of the processes of remodelling and mineralisation. If mineralisation kinetics were identical in trabecular and femoral cortical bone, our model predicts that the lower turnover rate in cortical bone would result in a broadening of the BMDD. This broadening is analogous to the broadening of the age distribution of humans when life expectancy increases: the human population becomes simply redistributed on a larger variety of possible ages. Our experimental observation that peak width in Fm.BMDD and Tb.BMDD is not significantly different seems to be observed in other situations of changed bone turnover analysed in terms of the BMDD:

The effect of different bisphosphonates on bone turnover (Ebetino et al., 2011) and the BMDD (Roschger et al., 2014) have been extensively studied over the last twenty years. Depending on their potency bisphosphonates can reduce the rate of bone remodelling substantially (Fuchs et al., 2011; Misof et al., 2013). After 10 years of alendronate treatment the mineral heterogeneity was not significantly different from the trabecular reference BMDD (Roschger et al., 2010). In a separate study, 5 years of risedronate treatment resulted in a CaWIDTH not significantly different from its baseline value (Zoehrer et al., 2006).

Hypoparathyroidism is a condition of known low bone turnover. In vitamin-D treated hypoparathyroid patients the activation frequency in trabecular bone of iliac crest biopsies was reduced by a factor 4.5 in these patients compared to healthy individuals (Langdahl et al., 1996). While the BMDD in the hypoparathyroid patients was significantly shifted towards higher values, the heterogeneity of mineralisation around the peak was not changed compared to the reference trabecular BMDD (Misof et al., 2016).

Our study also provides a new perspective on the definition of the trabecular reference BMDD of healthy adults (Roschger et al., 2003). The reference BMDD is defined based on the average of 52 individual BMDDs, since it was found that in healthy adults the trabecular BMDD does not depend on age, gender, ethnicity and skeletal site. The reference value for CaWIDTH was reported to be 3.35 wt% Ca with a standard deviation of ±0.34 wt% Ca (Table 2). Translating this standard deviation in a variability of turnover using Fig. 6 in Ref. (Ruffoni et al., 2007) shows that the corresponding variability in turnover would range from −40% to +10% of the average value, a remarkably narrow range having in mind the natural variability of turnover across bone sites and individuals, which is expected to be much larger (Parfitt, 2002). The fact that the trabecular BMDD is so consistent despite large variability of turnover rate suggests that another mechanism compensates the influence of turnover rate on mineral heterogeneity in trabecular bone.

The BMDD of cortical bone measured in transiliac bone samples was found to have a much lower CaPEAK value (Misof et al., 2014a) compared to the CaPEAK value of cortical bone at the midshaft femur (Fig. 1). Consequently, a reference cortical BMDD (independent of skeletal site) of healthy adults cannot be defined, in contrast with the reference trabecular BMDD. The similarity of CaWIDTH in Tb.BMDD and CaWIDTH in Fm.BMDD is inconsistent with what one would expect if differences in mineral density distributions would result from differences in turnover rate alone. If only turnover rate would differ, our model shows that Fm.BMDD would have a broader peak than Tb.BMDD, which is not what is observed.

Data from the control group of the risedronate study in Ref. (Zoehrer et al., 2006) shows that three years of treatment with calcium and vitamin D supplementation (with unknown consequences on the turnover rate) only resulted in a statistically significant shift of the BMDD peak by 0.73 wt% Ca (Borah et al., 2004; Zoehrer et al., 2006; Fratzl et al., 2007). This shift in peak location, presumably due to the extra calcium intake of the patients, was not associated with a change in BMDD peak width compared to baseline values.

A possible interpretation of our observations that the heterogeneity of mineral content around the mean seems rather uniform, is that the processes of remodelling and mineralisation are interlinked (see Eqs. (D.1), (D.2)). Indeed, our model shows that if mineralisation rate scales with turnover rate, then the degree of heterogeneity of mineral density around the mean, measured in particular by CaWIDTH, is independent of turnover rate.

A turnover-limited mineralisation kinetics suggests a limitation imposed by the transport and recycling of minerals. If mineral transport is spatially restricted, minerals that are embedded into bone matrix during bone mineralisation may come mostly from the recycling of minerals freed from the bone matrix during bone resorption (Roschger et al., 2020). The abundance of these recycled minerals depends on turnover rate.

In the studies (i)-(v) mentioned above the width of the BMDD peak was virtually unchanged, but the peak position typically shifted significantly. In our study an additional shift along the calcium content axis was needed to match the mineralisation kinetics of trabecular and femoral cortical bone. These shifts in the BMDD peak location suggest the existence of an additional calcium reservoir in bone that may be filled to different levels in trabecular and cortical bone during early stages of mineralisation. Indeed, the shift in calcium density between the mineralisation kinetic laws c(t) of trabecular bone and of femoral cortical bone is generated during the first 1–2 years of mineralisation, see Fig. 7. Other works have proposed the existence of mineral reservoirs in association with the lacuno-canalicular network of the osteocytes (Canè et al., 1982; Busse et al., 2010). In the context of a revived discussion about osteocytic osteolysis (Teti and Zallone, 2009; Bonewald, 2011), the osteocytes are thought to play an active role in the mineralisation process by modifying the perilacunar and pericanalicular matrix (Kerschnitzki et al., 2013; Hesse et al., 2015). In a recent study on human osteons, the local Ca content measured by qBEI was spatially correlated with the density of the lacuno-canalicular network measured by confocal microscopy. In line with the hypothesis of a mineral reservoir associated with the canalicular network, regions with a dense network were found to have a higher mineral content (Roschger et al., 2019).

Two alternative explanations for the difference between Tb.BMDD and Fm.BMDD should be mentioned. The first relies on the idea that the amount of calcium that can be incorporated into the bone matrix is limited and that this limitation may reduce CaWIDTH in low turnover scenarios (Misof et al., 2012, Misof et al., 2014a). Our previous computational analysis (Buenzli et al., 2018) exploring how a reduced capacity of mineral uptake influences the BMDD in late stages of mineralisation showed indeed that peak broadening is suppressed in low turnover rate scenarios. However, this computational analysis also predicted an increase in the BMDD's skewness which is not observed experimentally (Buenzli et al., 2018). The second possible explanation is based on the lack of reliable data about turnover rates depending on different skeletal sites and, in particular, different bone types (i.e., cortical vs. trabecular) (Parfitt, 2002). If femoral cortical bone and trabecular bone are turned over with more similar rates than what is usually assumed, then the similarity in the shape of the BMDDs of these different bones is also less surprising.

A strength of the quantitative hypothesis-testing approach used in this work is to be able to investigate the influence of various biological processes in isolation. The mathematical model we based our analysis on uses very general principles of mass conservation laws in which bone is created, matures, and is resorbed. The manner in which these processes occur was subject only to recovering the experimentally measured BMDDs. Nevertheless, the model makes two fundamental assumptions most important for the interpretation of our results. First, it is assumed that the temporal incorporation of minerals in bone matrix follows the same kinetics for bone formed at different location and different times in the same qBEI image. This assumption allows the definition of a mineralisation rate q(c), which can be different depending on bone type, but is otherwise the same for each newly formed bone packet in this bone type. While osteocytes may interfere with the time course of mineral accumulation during the mineralisation of bone matrix (Teti and Zallone, 2009; Atkins and Findlay, 2012; Belanger et al., 1967; Roschger et al., 2019), most measurements of mineralisation kinetics so far suggest that mineral accumulation remains a good proxy of tissue age within a same qBEI scan. Second, our analysis of adult BMDDs assumes that these mineral distributions are in steady state, i.e. they would not evolve in time under the given conditions of remodelling. This assumption seems justified considering that the age of the individuals is between 48 and 56 years old. However, this steady-state assumption may have to be revisited when interpreting BMDD data of individuals afflicted recently by certain bone disorders, or data from children due to their growth (Fratzl-Zelman et al., 2009b).

The mathematical model presented in this paper relies on the rigorous principles of the balance of elementary mineralised tissue volumes under remodelling and mineralisation processes. We applied these principles by considering a spatial scale that contains a similar level of detail as the experimental data. The model therefore includes mineral heterogeneity, but it does not include further variables such as spatial dependences. More detailed explorations of the influence of space and geometry (Buenzli et al., 2013; Lerebours et al., 2016; Alias and Buenzli, 2018) and the interplay with mechanics for the regulation of mineralised tissues (Huiskes et al., 2000; Gerhard et al., 2009; Dunlop et al., 2009; Christen et al., 2014) are of high interest for future works. The microarchitecture of trabecular bone gives rise to complex strain patterns, which could influence not only the remodelling, but also the mineralisation process. The consideration of bone tissues that retain traces of the history of their mode of growth in the distribution of mineral density is important to be able to analyse cortical BMDDs from the iliac crest, from children (Fratzl-Zelman et al., 2009b), and from diseased states.

In summary, our findings suggest that trabecular and femoral cortical bone are more different than usually assumed. The common view is that both bone types are made up of bone lamellae, but that there is a difference in the spatial arrangement of the lamellae, and a difference in the rate at which the bone is turned over. Our conclusion, drawn from quantitative hypothesis testing, that the mineralisation kinetics is different in these bone types suggests that this common view needs to be supplemented by the notion that minerals are incorporated in the two bone types differently. This conclusion is further corroborated by qBEI measurements of very low turnover rate trabecular bone, which exhibits a distinctively lower value of CaPEAK than that of femoral cortical bone despite having a similar low turnover rate (Misof et al., 2012, Misof et al., 2016). A possible reason for why trabecular bone and femoral cortical bone may exhibit differences in mineralisation kinetics is that these bone types have different mechanical requirements (Currey, 2006; Fratzl-Zelman et al., 2009a).

CRediT authorship contribution statement

Chloé Lerebours: Conceptualization, Methodology, Writing - review & editing, Writing - original draft, Software, Validation, Formal analysis, Investigation, Visualization. Richard Weinkamer: Conceptualization, Methodology, Formal analysis, Writing - review & editing, Supervision, Project administration, Funding acquisition. Andreas Roschger: Resources, Data curation, Methodology, Formal analysis, Writing - review & editing, Visualization. Pascal R. Buenzli: Conceptualization, Methodology, Formal analysis, Writing - review & editing, Supervision, Project administration, Funding acquisition.

Transparency document

Transparency document.

Declaration of competing interest

None.