Bone turnover and mineralisation kinetics control trabecular BMDD and apparent bone density: insights from a discrete statistical bone remodelling model

The mechanical quality of trabecular bone is influenced by its mineral content and spatial distribution, which is controlled by bone remodelling and mineralisation. Mineralisation kinetics occur in two phases: a fast primary mineralisation and a secondary mineralisation that can last from several months to years. Variations in bone turnover and mineralisation kinetics can be observed in the bone mineral density distribution (BMDD). Here, we propose a statistical spatio-temporal bone remodelling model to study the effects of bone turnover (associated with the activation frequency \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {Ac.f}$$\end{document}Ac.f) and mineralisation kinetics (associated with secondary mineralisation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_\textrm{sec}$$\end{document}Tsec) on BMDD. In this model, individual basic multicellular units (BMUs) are activated discretely on trabecular surfaces that undergo typical bone remodelling periods. Our results highlight that trabecular BMDD is strongly regulated by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {Ac.f}$$\end{document}Ac.f and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_\textrm{sec}$$\end{document}Tsec in a coupled way. Ca wt% increases with lower \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {Ac.f}$$\end{document}Ac.f and short \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_\textrm{sec}$$\end{document}Tsec. For example, a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {Ac.f}=$$\end{document}Ac.f= 4 BMU/year/mm3 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_\textrm{sec}$$\end{document}Tsec = 8 years result in a mean Ca wt% of 25, which is in accordance with Ca wt% values reported in quantitative backscattered electron imaging (qBEI) experiments. However, for lower \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {Ac.f}$$\end{document}Ac.f and shorter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_\textrm{sec}$$\end{document}Tsec (from 0.5 to 4 years) one obtains a high Ca wt% and a very narrow skew BMDD to the right. This close link between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {Ac.f}$$\end{document}Ac.f and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_\textrm{sec}$$\end{document}Tsec highlights the importance of considering both characteristics to draw meaningful conclusion about bone quality. Overall, this model represents a new approach to modelling healthy and diseased bone and can aid in developing deeper insights into disease states like osteoporosis.


Introduction
Bone is a living tissue that undergo continuous repair, renewal, and adaptation to its biochemical and mechanical environment.The bone remodelling process is essential for removing microcracks accumulated in the bone matrix due to dynamic loading and to regulate mineral homeostasis Edmund Pickering, Vittorio Sansalone, David Cooper and Peter Pivonka have contributed equally to this work.and hematopoiesis (Parfitt 1994).This adaptive behaviour is achieved through the bone remodelling process, which is performed by the basic multicellular units (BMUs) as introduced by Frost (1969).The bone remodelling process on trabecular bone is activated at the bone surface, and concerted action of osteoclastic bone resorption followed by osteoblastic bone formation controls bone turnover (Parfitt 1984(Parfitt , 2002)).During the bone mineralisation process the organic collagenous matrix deposited by osteoblasts becomes mineralised.The latter process is regulated by the mineralisation kinetics that exhibits two distinct phases: a fast primary mineralisation phase lasting for several days to a few weeks; and a secondary mineralisation phase that can last from several months to years (Bala et al. 2013).Given the nature of the mineralisation process, the mineralisation distribution depends strongly on the rate of bone turnover (Roschger et al. 2020), i.e. the number of BMUs being recruited, which can be described by the BMU activation frequency ( Ac.f ).While bone diseases can affect different levels of bone structure, such as trabecular architecture in osteoporosis (Wehrli et al. 2001;Parfitt et al. 1983) and degree of mineralisation in osteomalacia (Roschger et al. 1998), the inhomogeneous mineral content and its spatial distribution on a microscopic scale are major determinants of the mechanical quality of trabecular bone (Roschger et al. 2008).The spatial heterogeneous distribution of mineral determines the bone mineral density distribution (BMDD) (Boivin andMeunier 2002, 2003).Roschger et al. (2020) defined the BMDD as the distribution of calcium content and described it as the "fingerprint" of bone, since it can distinguish healthy from pathological bone tissue at the scale of the bone matrix (Ruffoni et al. 2007) (see Fig. 1), regardless of skeletal site, sex or ethnicity (Roschger et al. 2003).Parameters that are typically extracted from the BMDD distribution curves are (i) calcium peak ( Ca PEAK ), i.e. the value of calcium weight percertage at the peak of the BMDD curve, (ii) calcium width ( Ca WIDTH ), i.e. width of the BMDD frequency curve at half of the peak of the BMDD, (iii) mean calcium weight ( Ca MEAN ).
Several mathematical approaches aim to include the effects of bone turnover and bone matrix mineralisation into models of bone remodelling.In the following, we review the three major types of models and outline their pros and cons.
1. Continuous spatial averaging approach: In these models the bone remodelling process is described for a representative volume element (RVE), i.e. a typical trabecular bone RVE ≈ 5 mm 3 (Blöß and Welsch 2015).Usually, this approach formulates bone cell population models and the mineralisation process on this RVE while considering mechanobiological feedback (Martínez-Reina and Pivonka 2019; Martínez-Reina et al. 2021a).These types of models assume that several active BMUs in the RVE control bone turnover and density.However individual BMUs are not modelled explicitly, but only the accumulated action of osteoclasts and osteoblasts (of all active BMUs) control the changes in bone matrix volume fraction (or bone volume to total volume, known as BV/TV) and the average degree of mineralisation, i.e. ash fraction, in the RVE.These models have been typically applied in a pharmacokinetics-pharmacodynamics (PK-PD) modelling context to study the effects of different bone drug treatments on changes in BV/TV and mean mineralisation for patients with osteoporosis (Martínez-Reina et al. 2021a, b;Calvo-Gallego et al. 2022).We note that these models do not have the ability to recreate the BMDD and are not able to directly investigate the effects of changes in Ac.f and the spatial heterogeneity of the bone matrix.2. Continuous BMDD balance approach: Ruffoni et al.
(2007) described the changes in BMDD due to mineralisation as a "flow" from low to high values of the mineral content.They state that the area of the BMDD curve, i.e. calcium percentage by weight (Ca wt%) versus frequency in a constant bone volume, remains unchanged during the remodelling process.In fact, they conclude that the effect of remodelling on the BMDD is seen by a flow from lower Ca wt% values towards higher Ca wt% values (osteoblast action), while some bone volume is lost due to a "leakage" of the flow (osteoclast action).These effects were described using a BMDD balance equation in analogy with mass balance that can be mathematically expressed as a reaction-advection equation.Using this formulation, they investigated different mineralisation laws and their effects on BMDD in the RVE (Ruffoni et al. 2007).In subsequent work, Buenzli et al. (2018) extended this model to limit mineralisation to various maximum calcium capacities of bone.Results of this extended model show that an abrupt stopping of mineralisation near a maximum calcium capacity induces a pile-up of minerals in the BMDD statistics that is not observed experimentally.Using a smooth decrease in mineralisation rate, imposing low maximum calcium capacities helps to match peak location and width of simulated low-turnover BMDDs with experimental BMDDs.However, this study resulted in a distinctive asymmetric peak shape of the BMDDs (Buenzli et al. 2018).We note that even though these models are able to investigate the BMDD, no conclusions can be made regarding the bone microstructure.3. Discrete statistical approaches: Martin (1984Martin ( , 1991) ) showed that calculation of porosity and bone volume fraction are the same for two-dimensional (2D) and three-dimensional (3D) stereological approaches.In particular, discrete models allow us to calculate the porosity evolution, assuming perfect cylindrical osteonal remodelling perpendicular to the cortical cross section.Using discrete models of cortical bone remodelling Martin (1984Martin ( , 1991) ) investigated various patterns of osteonal remodelling including overlap feedback, i.e. new secondary osteons have the tendency to remodel existing osteons, and random remodelling on the evolution of cortical porosity (Hazelwood et al. 2001;Nyman et al. 2004).Moreover, Heaney (1994); Heaney et al. (1997) used a discrete approach to model bone remodelling by assuming a certain number of BMUs in the RVE which could be integrated over a remodelling period to study the transient behaviour of bone remodelling.Later, Thomsen et al. (1994) developed a stochastic model of the remodelling process for human vertebral trabecular bone.Their simulations predicted the long-term effects of changes in the remodelling process on bone mass, trabecular thickness, and perforations of trabeculae.However, previous models do not include the mineralisation process and therefore cannot assess the calcium distribution.
Overall, even though the BMDD has great value to assess bone quality and identify bone diseases, the available bone remodelling models are not capable of recreating the spatial calcium heterogeneity and investigate the effects on BMDD and apparent density of factors driving bone remodelling.Therefore, we propose a discrete statistical two-dimensional computational model of the bone remodelling process.
Based on the works of Martin (1984Martin ( , 1991)); Heaney (1994), this remodelling model of the trabecular bone is discrete both in space and time, but differently than previous discrete models, we include every phase of the remodelling process, from resorption to mineralisation.In this work, we focus on studying the effects of BMU Ac.f and mineralisation kinet- ics on the BMDD, as well as the bone microstructure.We apply this new approach to healthy and diseased bone, i.e. osteoporosis, and give insights on the contribution of bone remodelling factors on the overall bone quality.

Discrete remodelling model for trabecular bone
All parameters used in the model can be found in Table 1, where symbols are based on the standardised nomenclature and symbols for bone histomorphometry (Dempster et al. 2013).

Structure and composition of bone tissue
Let us consider a 2D region Ω RVE of trabecular bone, small enough with respect to the whole trabecular compartment of the bone it belongs to.This region describes a slice of a Representative Volume Element (RVE) (around 5 mm 3 in trabecular bone (Cowin 2001)).Hereinafter, following a stereological assumption, we will focus on a 2D description of bone and consider that our analysis is consistent with a 3D description.Therefore, bone composition will be described in terms of areas and area fractions of its constituents, that are assumed to be representative of their 3D counterparts-volumes and volume fractions, respectively.At the tissue scale, trabecular tissue is made up of bone matrix and marrow, occupying the regions Ω Bm and Ω Ma , respectively.It follows that: where A RVE , A Bm , and A Ma are the areas of the whole tra- becular tissue, and of the bone matrix and marrow regions, respectively; moreover, t denotes the time.
Equation (1) can be normalised with respect to A RVE , and written in terms of area fractions: where f Bm and f Ma are the area fractions of bone matrix and marrow space, respectively, at the tissue scale.
Let p be the 2D vector that denotes the position of a point in Ω RVE .Since our model is discrete in space and time, every point p actually refers to a square pixel with side length equal to , here called the model resolution in space, and measured in μ m.All p ∈ Ω Bm is considered to be part of the bone matrix, which basically consists of organic matter (mostly collagen), mineral, and water.Thus, neglecting the other minor (1) constituents of bone matrix, the composition of each point of the bone matrix can be written as: where o , m , and w are the area fractions of organic matter, mineral, and water, respectively, at the pixel scale.

Mass densities, ash factor, and calcium content
Let be the mass density of the bone matrix.It can be expressed as: where o , m , and w stand for the mass densities of organic matter, mineral, and water, respectively.All these mass densities are assumed not to change over time.
For the sake of comparison with experimental observations, it is useful to compute the apparent mass density of bone matrix: where Nb RVE is the number of points p ∈ Ω RVE .
The mineral content in the bone matrix is usually characterised by the ash fraction , which describes the ratio between the mineral mass m m and dry mass m dry (the sum of mineral and organic mass), reading: Provided the ash factor, computing the calcium content is straightforward.The calcium content is defined as the mass of calcium m Ca per dry mass, reading: . where we dropped the dependency of the variables on space and time for the sake of clarity.The coefficient on the righthand side of the above equation is obtained through the stoichiometric formula for hydroxyapatite.Hereinafter, we will denote Ca MEAN and Ca PEAK the mean and peak value of the calcium content distribution, respectively (Roschger et al. 2003) (see Fig. 1).

Remodelling process
Let us assume that the 2D region Ω RVE is experiencing remodelling and let B.Pm describe the bone matrix perim- eter, i.e. the set of points belonging to the perimeter of Ω Bm , which is analogous to the trabecular surface in a 3D geometry.Remodelling is performed by basic multicellular units (BMUs), that are activated on B.Pm with a frequency Ac.f .A remodelling event is depicted in Fig. 2. Remodel- ling starts with BMU activation at the point p S (Fig. 2A) and proceeds through the resorption (Fig. 2B), reversal (not shown), formation, and mineralisation (Fig. 2C) phases.First, we assume that the osteoclasts dig a semi-circular area of bone for a resorption period Rs.P , as proposed by Parfitt (1994) in 1994.This is followed by a reversal period Rv.P in which the transition from resorption to formation takes place.Subsequently, we consider that osteoblasts lay down osteoid in semi-circular layers for a period FP .Finally, mineralisation takes place as new bone is laid down, after a mineralisation lag time Mlt .The mineral content increases during mineralisation, as mineral ions migrate from the interstitial fluid into the bone matrix to form hydroxyapatite ( Ca 10 (PO 4 ) 6 (OH) 2 ).These crystals nucleate within and outside the collagen fibrils, displacing the water present in the bone matrix (Bala et al. 2013).
The time course of the radius of the semi-circular region where remodelling takes place is depicted in Fig. 3.Each BMU is active during a time interval ( 7) Ac.P = Rs.P + Rv.P + FP = T 3 − To .For each remodel- ling event, two local time counters are introduced: b to follow the evolution of the active BMU across the activation, resorption, reversal, and formation phases; and m to follow the mineralisation phase.We assume that osteoclastic bone resorption and osteoblastic bone formation are linear processes.The dashed line refers to unbalanced remodelling, namely when bone resorption prevails on bone formation-see Sec.2.5.4 for more details.
It is important to note that we do not model the cellular dynamics related to bone remodelling.We assume that individual BMUs within the RVE have reached a steady-state configuration, i.e. all bone cell populations have reached their steady-state distributions (Buenzli et al. 2014).

Remodelling algorithm
The iterative algorithm of our discrete model is shown in Fig. 4. For every simulation, we input the trabecular bone geometry from micro-CT images; the initial Ca wt% in the Ω RVE ; the hemiosteonal geometry as described in Fig. 3; the mineralisation law; and the activation frequency Ac.f (see Table 2 for more details).From the Ac.f and total simula- tion time T, we identify the total amount of BMUs to be created.We initiate the simulation at time t = 0 .When t is equal to n times T rep (counter of remodelling events times the periodicity of BMU activation events) a new BMU is created and we initiate a local time counter b that follows the evolution of the active BMU.While the BMU is active, it goes though every stage of the bone remodelling process: resorption, reversal, formation and mineralisation.At every time b the output variables var are updated, where var includes area fractions and the calcium concentration ( var = f Bm , f Ma , o , w , m , Ca wt% ).It is important to note that multiple BMUs can be active simultaneously, since each BMU holds their own local time counter b .Once a point p goes through formation it can be mineralised, and a mineralisation local time counter m is initiated for every single formed point p.The mineralisation level evolves through a mineralisation law, described in Sec.2.5.5.Points in Ω Bm that do not belong to any active BMU and still have not attained the mineralisation threshold have their mineralisation level and mineralisation local time counter m updated at the end of every iteration.Once all active BMUs and all bone have gone through their remodelling phase the time t is updated, and the process is repeated until t reaches the total simulation time T.

Remodelling phases
The different phases of the remodelling process and their implementation in our model are described in this section.

Activation phase
The factors that initiate bone remodelling process are not completely understood, but BMU activation is believed to occur partly due to biomechanical demands, e.g.removal of microcracks from the bone matrix-hereinafter referred to as targeted remodelling-and due to biological demands, e.g.calcium and phosphate homeostasis-referred to as random remodelling (Wss 2001).In our model, the activation phase is governed by the activation frequency Ac.f , i.e. the number of BMUs activated per unit time and unit bone area, that describes the intensity of remodelling.Ac.f is meant to be an input parameter of our model and is measured in BMU/ mm 2 /year.The total number of BMUs activated during a time interval Δt in the bone region Ω RVE is given by: (8) N = Ac.fΔt A RVE .Fig. 4 Discrete model iterative algorithm.After initiation of parameters and geometry, the simulation starts at time t = 0 .If t is equal to n times T rep (counter of remodelling events times the periodicity of BMU activation events) a new BMU is created and a local time counter b starts to follow the evolution of the active BMU.Every BMU goes through every stage of the remodelling process.The duration of each event is given by the local time points ( 1 , 2 , 2 ), as explained in Fig. 3. Once a point p goes through formation it can be mineralised, and we initiate a mineralisation local time counter m for every single formed point p. p ∈ Ω Bm that do not belong to any active BMU and still have not attained the mineralisation threshold have their mineralisation level and mineralisation local time counter m updated at the end of every iteration.Once all active BMUs and all bone have gone through their remodelling phase the time t is updated, and the process is repeated until t reaches the total simulation time T. The model output var includes f Bm , f Ma , o , w , m , and Ca wt% Our model can describe both random and targeted remodelling processes.In random remodelling, BMUs are randomly activated at any point of the boundary B.Pm .In targeted remodelling, BMU activation is related to the mineralisation of the bone matrix, as microcracks are more likely to appear at bone sites which are highly mineralised.Hereinafter, we will focus on targeted remodelling.
Let n be a counter of remodelling events.We assume that remodelling events are evenly distributed in time, the periodicity being T rep = 1∕(Ac.fA RVE ) .Thus, the n-th remodel- ling event starts at t n = n T rep .In order to activate the BMU, we search the Ω RVE for the regions on the boundary of the bone matrix featuring the highest average degree of mineralisation.More precisely, for each p ∈ B.Pm we define a search region around p as shown in Fig. 2A (the size of the search region is a model parameter) and compute the average value of within it, to be called ᾱ .The activation point of the n-th BMU, referred to as p n S , corresponds to the region featuring the highest value of ᾱ .If several regions present the same value of ᾱ , p n S is chosen randomly among those regions.
Eventually, the time counter n b is started and remodelling moves to next phase.

Resorption phase
Trabecular resorption may produce different patterns, e.g.trenches and pits.For the sake of simplicity, we assume resorption to produce hemiosteonal patterns (Parfitt 1994), i.e. semi-circular cavities with (hemiosteonal) radius R On -see Fig. 2B.
The resorption phase of the n-th BMU takes place over the local time interval [0.. 1 ] , corresponding to the global time interval [T n 0 ..T n 1 ] .Bone resorption starts at p n S (for n b = 0 ) and proceeds radially up to a distance R On from p n S (for n b = Rs.P ).Let r Rs be the current resorption radius, defining the limits of the current resorption region.We assume r Rs to evolve linearly in time, namely: The current resorption region is:

Reversal phase
The reversal phase links bone resorption and bone formation.It sets the stage for the subsequent phases of bone remodelling, allowing for the deposition of new bone tissue.
The coupling mechanisms between osteoclasts and osteoblasts responsible for generating an osteogenic environment remain poorly understood (Delaisse 2014).In our model, the reversal phase is a quiescent phase.
The reversal phase of the n-th BMU takes place over the . Throughout this phase, the radius of the remodelling region is equal to R On .

Formation phase
In the formation phase, osteoblasts lay down layers of unmineralised bone matrix called osteoid.Osteoid is composed of collagen and other proteins that provide the framework for mineralisation.Similarly to the resorption, in our model bone formation is assumed to proceeds radially from the hemiosteonal wall (i.e. the cement line) towards p n S -see Fig. 2C.
The formation phase of the n-th BMU takes place over the local time interval . )and proceeds radially towards p n S leading to complete ( n b = 3 ) or incomplete (  n b <  3 ) refilling of the resorption cavity.Let r F be the current formation radius, defining the limits of the current formation region.We assume r F to evolve linearly in time, namely: The current formation region is: Just after being laid down, new bone is completely unmineralised ( m = 0 ).This new bone is basically a hydrogel made of collagen and water.The collagen content of bone can vary depending on factors such as species, age, diseases, anatomical location, and bone type (cortical or trabecular).However, variations are generally limited.The volume fraction of the organic matter is assumed to be the same for all the bone matrix points and not to change over time.The value o = 1 3 has been used hereinafter (Vuong and Hellmich 2011;Martin 1984).Thus, in view of Eq. ( 3), the initial composition of each new point p ∈ Ω n F is set as: Bone resorption and formation are balanced in a healthy remodelling process.By contrast, reduced bone formation can be associated with pathological bone remodelling.In our model, healthy remodelling corresponds to the case where the hemiosteonal cavity-dug during the resorption phase-is entirely filled during the formation phase.Accordingly, unbalanced remodelling can be simulated by considering that the hemiosteonal cavity is only partially filled during the formation phase.The formation radius at the end of the formation phase can be expressed as u f R On , with u f being the underfilling coefficient.The latter is a model parameter taking values between 0 (complete refilling of the hemiosteonal cavity) and 1 (no filling).We implemented partial refilling in our model by stopping the formation phase as soon as r F (as per Eq. ( 11)) drops below u f R On .

Mineralisation phase
Every point of bone matrix goes through the mineralisation process.Mineralisation of osteoid (laid down in the formation phase) starts after a mineralisation lag time Mlt , during which collagen fibres within the osteoid organise so as bone can form properly.In our model, after a new point p ∈ Ω n F is laid down at (global) time t F (p) , an associated time counter m is introduced to follow its mineralisation, and it is started after Mlt from t F (p) , i.e. m (p, t) Once mineralisation has started in the osteoid, a continuous increase in its mineral content occurs.Mineralisation is described by the evolution of the ash factor according to a mineralisation law M( m ) , namely: (p, t) = M m (p, t) .As discussed in the introduction, primary mineralisation constitutes the first rapid phase of the mineralisation process.In this period, the ash factor is reported to reach the value = 0.45 (Hernandez et al. 2000).The mineralisation process is then followed by a slow and gradual increase in the mineral content over a longer period of time, i.e. the secondary mineralisation, which evolves to a maximum Ca concentration of 30 wt% (Currey 2004;Ruffoni et al. 2007).We tested all the mineralisation laws proposed by Ruffoni et al. (2007) to represent both stages of mineralisation.The most promising BMDDs were obtained through the double exponential and hyperbolic mineralisation laws.For the sake of simplicity, we describe the mineralisation law through a double exponential curve: where c 1 and c max are calcium contents, and t 1 and t 2 are char- acteristic times-see App. 1 for more details.It is important to note that the total time of the secondary mineralisation is defined by T sec , and the slope of the secondary mineralisa- tion curve, i.e. second exponential, is given by t 1 .
As mentioned in Sec.2.5.4, the area fraction of the organic matter is assumed to be the same for all points and not to change over time.Thus, mineralisation of a bone matrix point emerges as a perfectly balanced increase of its mineral area fraction m and decrease in water area fraction w .Using Eq. ( 6), the mineral area fraction can be written for each point p ∈ Ω Bm as: Eventually, since o is known and m (p, t) is provided by Eq. (15), Eq. (3) allows computing w (p, t).

Model implementation
The presented model was implemented and evaluated using MATLAB R2022a with suitable initial conditions to investigate changes in BMDD distribution due to variation in activation frequency and mineralisation kinetics.

Geometry
We define the Ω RVE geometry from and in-silico generated trabecular bone as typically observed by micro-CT imaging.The trabecular bone generated image has a resolution of 1μm × 1μm and trabecular thickness ranging from 200μ m ( 14) to 400μ m (Liu et al. 2010).We use an edge detection tool implemented in MATLAB R2022a to identify the bone perimeter B.Pm .The Ω RVE is binarised into bone and mar- row areas, and lacuno-canalicular porosity in bone matrix is neglected, which contributes about 3-5 % porosity (Cowin 1999;Ashique et al. 2017)

Initial state and model parameters
At the initial state, we assumed the ash factor uniform throughout Ω Bm and equal to 70%.Model parameters related to the remodelling process are adapted from the literature and reported in Table 2.Note that the resorption radius R On is defined randomly at each new BMU initiation according to a normal distribution with mean value ROn and standard deviation On .The total time of secondary mineralisation ( T sec ) and the activation frequency ( Ac.f ) are free param- eters (f.p.) of our model and will be discussed afterwards.Eventually, the time step in all simulations is equal to 1 day.

Secondary mineralisation characteristic time
The majority of works agree that the primary mineralisation takes about 5-10 days (Boivin et al. 2009), but discrepancies are found regarding secondary mineralisation time T sec .While some experimental studies state that secondary mineralisation takes place up to several months (Berli et al. 2017), several numerical studies use values of up to 10 years for secondary mineralisation (Martínez-Reina et al. 2008;Bala et al. 2007).Aiming analyse the mineralisation kinetics, we studied different values for T sec , as seen in left panel of Fig. 5.More details about the secondary mineralisation time T sec and the parameters of the mineralisation law M in Eq. ( 14) can be found in Appendix 1.
Even though the mineralisation laws are written in terms of the ash factor , the evolution of mineral fraction over time can be obtained by using Eq. ( 15).The m in function of is shown in right panel of Fig. 5.

Activation frequency
As T sec , Ac.f values vary greatly in the literature.On the one hand, Nyman et al considered Ac.f = 4 BMU/mm 2 /year for bone remodelling in premenopausal trabecular bone.On the other hand, the experimental study performed by Parfitt (1983) suggests that Ac.f in trabecular bone is equal to 18 BMU/mm 2 /year.Moreover, values of Ac.f also vary for path- ological cases.Indeed, the unbalance of bone resorption and formation seen in postmenopausal osteoporosis is said to be explained by a rise in Ac.f (Sambrook and Cooper 2006).

Results
Given the random nature of the proposed remodelling algorithm, we run fifteen simulations using the same initial geometry, as shown on the left panel of Fig. 9.For each simulation, the random generator seed was initialised based on the current time, resulting in a different sequence of random numbers.Unless otherwise stated, the results are presented in terms of the mean of all the simulations and the parameters presented in Table 2 have been used in the simulations.Ac.f  and T sec First, we study the effect of Ac.f and mineralisation kinetics in the healthy remodelling process.Figure 6 presents the average calcium content Ca MEAN and peak calcium content Ca PEAK , at T = 5 years for varying activation frequency Ac.f (ranging from 1 to 8 BMU/mm 2 /year) and secondary mineralisation periods T sec (ranging from 0.5 to 12 years).It can be noticed that Ca MEAN decreases as either of Ac.f or T sec increases.The standard deviation of Ca MEAN varies from 0.0012 to 0.0477 wt%, and is therefore not shown.Moreover, the Ca PEAK only decreases for high values of Ac.f and slow secondary mineralisation.The effects of Ac.f and T sec on the distribution of Ca wt% can be appreciated looking at the variations of the BMDDs with respect either Ac.f (Fig. 7) or T sec (Fig. 8).On the one hand, for a given mineralisation kinetics ( T sec = 8 years), the BMDD follows a negatively skewed distribution for lower values of Ac.f , and the peak of Ca wt% is about 30% (Fig. 7A-D).On the other hand, skewness of the BMDD decreases with higher values of Ac.f , and the peak of Ca wt% decreases towards 25% (Fig. 7E-H).

BMDD mean decreases with increase in
A similar pattern is seen for a constant Ac.f and varying mineralisation kinetics (Fig. 8).Assuming Ac.f of 4 BMU/ mm 2 /year, the BMDD follows a negatively skewed distribution for fast mineralisation kinetics, and the peak of Ca wt% is about 30% (Fig. 8A-D).Skewness of the BMDD Fig. 6 Mean and peak of calcium concentration ( Ca MEAN and Ca PEAK ) versus Ac.f for different mineralisation kinetics ( T sec ), evaluated at simulation time T=5 years Fig. 7 BMDD for mineralisation kinetics T sec = 8 years and Ac.f equal to A 1 BMU/mm 2 /year, B 2 BMU/mm 2 /year, C 3 BMU/mm 2 /year, D 4 BMU/mm 2 /year, E 5 BMU/mm 2 /year, F 6 BMU/mm 2 /year, G 7 BMU/mm 2 /year, and H 8 BMU/mm 2 /year decreases with longer times of second mineralisation, with the peak of Ca wt% decreasing towards 26% (see Fig. 8E-F).
As mentioned earlier, values of Ac.f and T sec vary greatly in the literature.For example, the experimental study performed by Parfitt (1983) suggests that Ac.f in trabecular bone is equal to 18 BMU/mm 2 /year.We performed five simulations using the suggested parameters ( Ac.f = 18 BMU/mm 2 /year and R On = 40μm).The initial geometry of the simulation is shown in the left panel of Fig. 9; the solutions obtained at T = 8 years using typical values used in our analysis and those suggested by Parfitt are shown in the middle and right panels of Fig. 9, respectively.The latter simulation shows large areas of interstitial bone that are never remodelled.A deeper insight on the BMDD obtained using the new parameters is shown in Fig. 10. Figure 10A shows the variation of the mean calcium content Ca MEAN with respect to the secondary mineralisation time T sec .Val- ues of Ca MEAN between 22% and 24% have been obtained for slower mineralisation kinetics, which are in accordance with Ca MEAN values seen in the literature (Lerebours et al. 2020).BMDD plots obtained for T sec = 6 years and 8 years are shown in Fig. 10B and C, respectively.BMDDs are negatively skewed and the skewness decreases for slower mineralisation kinetics.Moreover, the frequency peaks at the right of the BMDDs (Ca wt% = 30%) represent bone that is never remodelled in the presented simulations.

Bone density as function of Ac.f and T sec
We also analysed the variations of the apparent density app with respect to the activation frequency ( Ac.f ranging from 1 to 8 BMU/mm 2 /year) and mineralisation kinetics ( T sec rang- ing from 0.5 to 12 years)-see Fig. 11.As for Ca MEAN , results in Fig. 11 (left panel) show that app decreases as either Ac.f or T sec increases.On the other way round, low bone turno- ver combined with faster secondary mineralisation results in higher densities.
As shown in Fig. 11 (right panel), for a given Ac.f , there is a mineralisation kinetics that maintain a constant app through the remodelling process.For example, a combination of Ac.f = 4 BMU/mm 2 /year and T sec = 4 years results in small variations of app that remains about 0.68 g/cm 3 .

Bone remodelling and osteoporosis
Finally, we simulate different types of primary osteoporosis, such as post-menopausal (type I) and senile (type II) osteoporosis.From a biological point a view, post-menopausal   (Martin et al. 1998) and Ac.f = 6 BMU/mm 2 /year and (right) R On = 40μ m and Ac.f = 18 BMU/mm 2 /year as proposed by (Parfitt 1983).Both cases consider T sec = 8 years osteoporosis is associated with an unbalanced activity of osteoclasts and osteoblasts, namely when bone resorption prevails on bone formation.This imbalance can be magnified by a rise in Ac.f (Sambrook and Cooper 2006).Therefore, we simulate type I osteoporosis following the hypothesis that it is related to longer formation periods FP (Villanueva et al. 1966) or to high Ac.f (Sambrook and Cooper 2006) (Fig. 12A and B, respectively).On the other hand, senile osteoporosis is mainly characterised by a shift from osteoblastogenesis to predominant adipogenesis in the bone marrow that results in a reduction in bone formation (Duque and Troen 2008).Accordingly, we describe senile osteoporosis following the hypothesis that it occurs through an underfilling of semiosteons or that no formation takes place, which is usually named as uncoupling (Delaisse et al. 2020) (Fig. 12C and D, respectively).
Since osteoporosis is characterised by a low bone mass and micro-architectural deterioration associated with a negative net balance between bone formation and resorption during bone remodelling (Langdahl et al. 2016), we analysed the evolution of the bone matrix ( f Bm ) and marrow ( f Ma ) fractions, and of the apparent density ( app ).We first performed a simulation in physiological conditions ( Ac.f = 4 BMU/mm 2 /year and T sec = 4 years).The configuration obtained after 5 years of physiological remodelling was used as initial configuration for the simulations describing pathological remodelling.We simulated different scenarios leading to osteoporotic bone, by modifying the parameters used for physiological remodelling: the High FP scenario corresponds to FP = 475 days (Villanueva et al. 1966); the High Ac.f scenario corresponds to Ac.f = 8 BMU/mm 2 /year; the Underfilling scenario corresponds to u f = 50 % ; and the No filling scenario corresponds to u f = 1 .High FP and High Ac.f scenarios refer to postmenopausal osteoporosis.Sce- narios Underfilling and No filling scenarios refer to senile osteoporosis.All these simulations describe the evolution of When the formation phase is prevented and BMUs are not filled (No filling scenario), a drastic loss of bone mass and an important decrease in app can be observed.This is a severe case of osteoporosis where the microstrucure is severely affected (Fig. 12D), and the bone recovery is extremely unlikely, since it may include perforation of trabeculae.For the other cases, variations in bone matrix fraction and apparent density are less important, but still noticeable when compared to the control.The Underfilling scenario shows a marked decrease of both bone matrix fraction f Bm and apparent density app .Both High Ac.f and High FP scenarios show a limited decrease of f Bm .However, the latter show a pronounced decrease of app , comparable to that of the Underfilling scenario.
When looking at the BMDD, it can be observed that the Ca MEAN decreases for the postmenopausal case with high Ac.f (see Fig. 14 left).This decrease is expected, since an increase in Ac.f results in a less mineralised tissue.Our numerical results agree with experimental observations that show a shift of the BMDD towards lower values of Ca content in postmenopausal cases (Zoehrer et al. 2006) (see also Fig. 1).Furthermore, for senile osteoporosis, the Ca concentration has a tendency to increase, since little to no bone is being formed while bone continues to mineralise (Fig. 14 right).

Discussion
The current model considers the process of mineralisation of the bone matrix, differently than existing discrete remodelling models, such as the ones proposed by Martin (1984Martin ( , 1991)); Nyman et al. (2004).As discussed by many authors, computational models that neglect bone matrix mineralisation are not able to predict increases in apparent bone density (or intrinsic bone density) for cases when bone remodelling is significantly reduced, such as for example in the use of anti-resorptive drugs (Scheiner et al. 2014;Martínez-Reina and Pivonka 2019).Values used for the activation frequency Ac.f in our numerical simulations are lower than values suggested by some experimental studies (Parfitt 1983).These lower values could be because the width of hemiosteons may vary with the geometry of trabeculae.If hemiosteons are smaller, a higher Ac.f is needed to achieve the same bone turnover.The value of 18 BMU/mm 2 /year was proposed by Parfitt (1983) for hemiosteons having a radius of 40μ m.In our simula- tions, we used a hemiosteon radius R On = 100μ m.Therefore, Ac.f should be considerably smaller, as used in the current simulations, which agrees with the values of Ac.f used by Nyman et al. (2004).
We found that BMDD depends both on Ac.f and T sec shown in Figs. 7, 8.This finding supports the conceptual model of Boivin andMeunier (2002, 2003) which suggests these two factors to affect the BMDD.Similar trends in BMDD were found by Ruffoni et al. (2007) et al using the continuous BMDD balance approach.However, unlike in Ruffoni et al. (2007), the BMDD histograms of our simulations were not sensitive to the type of mineralisation curve (i.e.hyperbolic versus exponential) as long as T sec and the maximum level of mineralisation were kept similar.
For situations where secondary mineralisation is fast (i.e.T sec equal to one year or below), most of bone is highly mineralised, and this leads to a negatively skewed BMDD histogram (Fig. 8A-B).Therefore, although Martin (1991) states that variation in Ac.f alone can be responsible for bone remodelling an define bone gain and loss, the mineralisation kinetics play an essential role in determining the BMDD.In accordance with Lerebours et al. (2020), we found that Ac.f and T sec have a coupled effect on the BMDD.
For lower activation frequency or fast secondary mineralisation, the BMDDs are strongly negatively skewed.Skewness of BMDDs decreases for higher activation frequency or slow secondary mineralisation, as observed through quantitative backscattered electron imaging (qBEI) (Roschger et al. 1998).qBEI experimental results on trabecular bone suggest values of Ca MEAN ranging between 22 and 24 wt% (Lerebours et al. 2020).Our numerical results show a Ca MEAN slightly higher (between 24 and 26 %).This could be explained by some approximations of the model.Firstly, only a two-dimensional slice of trabecular bone was considered, which leads to some inner trabecular regions having never been remodelled, because BMUs are activated on bone surfaces only.Secondly, BMU activation in the bulk of the trabeculae are not considered by our model.However, bulk remodelling within trabeculae-creating cylindrical osteons similar to those observed in cortical bone-does occur (Sato et al. 1986), and could serve as a means of targeting areas of higher Ca content within thicker trabeculae.
Unlike Ca MEAN , the apparent bone density only changes slightly for different Ac.f and T sec .In the studies cases, app varies from 0.6 to 0.75 g/cm 3 , which agrees with density values of trabecular bone seen in the literature (Zioupos et al. 2008;Öhman-Mägi et al. 2021) The osteoporotic remodelling patterns highlight that even though the apparent density for type I and type II osteoporosis is comparable, the former is due to a high percentage of newly formed bone, while the latter is due to general bone loss.Therefore, our results suggest that analysing only the apparent bone density is not be sufficient to understand the nature of bone pathology.Moreover, besides presenting a lower app , the case with higher Ac.f also shows a shift in Ca MEAN , as seen experimentally (Zoehrer et al. 2006).These results suggest that type I osteoporosis is mainly due to high Ac.f rather than longer FP.
Finally, our study has a number of limitations.Firstly, no mechanical loading has been included in our simulations.However, mechanical loading has been shown to potentially affect mineralisation kinetics (Tourolle né Betts et al. 2020).Also, mechanical loading may induce microcracks in the bone matrix which could serve as a target for BMU initiation.Here, we used higher mineralised bone tissue as a target for BMU activation.Secondly, BMU resorption, reversal and formation events were assumed to vary linearly, and were implemented in an ad hoc fashion with no biological mechanism considered.Lastly, trabecular BMUs were described by Parfitt (1994) as spatio-temporal systems where osteoblasts follow osteoclasts as they progress along the bone surface.This was depicted as equivalent to half a cortical BMU  2020) depict a trabecular BMU as progressing and creating a "resorption track".Our modelling is not inconsistent with trench-like BMUs progressing perpendicularly to the 2D cross section of trabecular bone considered in our study, leading to semicircular patterns in the cross section.However, it cannot describe BMUs progressing on in other directions, namely along the bone surface in the plane of the cross section.With respect to the cross-sectional shape of trabecular BMUs, a recent histological 2D study of osteopontin-stained cement lines in human vertebral samples, revealed a lens-shaped morphology, deepest in the centre and progressively thinner towards their peripheries (Lamarche et al. 2022).While trabecular remodelling has been extensively studied, there is relatively little direct 3D morphological BMU data published.Serial block face imaging revealed trabecular remodelling events characterised by indentations on trabecular surfaces, which is relatively similar to the idealised forms we modelled (Slyfield et al. 2012;Tkachenko et al. 2009).However, more extended trench-like forms have also been reported (Tkachenko et al. 2009).Kragstrup and Melsen (1983) described trabecular BMUs from human samples as "broad bands" or "plate shaped" with extended trench-like structures, rather than punctate forms.Furthermore they observed that hemiosteons could be curved, branched and even cylindrical, which is not taken into account in our work and should be studied in future works.
The ultimate aim of our modelling approach is to also add mechanics to this model.For instance, introducing information about dynamics of loading cycles and material behaviour may allow to investigate the coupling between BMDD and stress patterns in the bone.In turn, this will allow to analyse and track accumulation of microcracks in the bone matrix, and ultimately estimate risk of bone fracture, which would depend on BMDD.

Conclusions
In this work, we propose a discrete statistical model of BMU remodelling in trabecular bone.Besides its geometrical algorithmic nature, this discrete model allows the study of the evolution of bone microstructure, apparent bone density, and bone mineral density distribution as a function of BMU activation frequency and time of secondary mineralisation.The major findings of our numerical simulations are: • BMDD is strongly influenced both by activation frequency of BMUs and mineralisation kinetics.• Normal BMDD histograms for healthy equilibrium bone remodelling situations are obtained for different combinations of Ac.f , and T sec values.
• For equilibrium cases Ca MEAN depends on the value of Ac.f with higher Ac.f giving lower Ca MEAN values.• Apparent bone density and BMDD should be mutually analysed to understand the nature of bone diseases, especially osteoporosis.
Future model developments will include mechanical aspects of bone remodelling and mechanobiological regulation of BMUs.

Appendix A Mineralisation law parameters
Table 3 presents the parameters for the double-exponential mineralisation law described in Eq. ( 14) for different secondary mineralisation kinetics.Note that the only varying parameter is the characteristic time t 1 , which is mainly responsible for the "slope" of the second exponential curve, i.e. the secondary mineralisation.By contrast, this parameter does not affect noticeably the first

Fig. 1
Fig. 1 Comparison between BMDDs of trabecular bone from a healthy population (data from Roschger et al. (2008)), from a postmenopausal osteoporotic (PmOP) population, and from an osteomalacia population (data fromZoehrer et al. (2006)).Ca PEAK indicates the peak in the distribution, Ca MEAN the average degree of mineralisation, and Ca WIDTH the width at half of the peak of the BMDD

Fig. 2
Fig. 2 Diagram showing the development of a BMU on a cross section of trabecular bone A RVE (white and red represent bone and marrow, respectively): A Activation process giving position of a new

Fig. 3
Fig. 3 Hemiosteonal geometry and time course of the radius of the semi-circular region where remodelling takes place.Besides the global time axis t, local time axes b and m are introduced to follow the progression of the BMU and mineralisation, respectively.BMU is activated and resorption phase starts at local time b = 0 , correspond- [ 2 .. 3 ] , corresponding to the global time interval [T n 2 ..T n 3 ] .Bone formation starts at the cement (

Fig. 5
Fig. 5 Left panel: Effect of secondary mineralisation characteristic time ( T sec ) on the ash factor ( ) and calcium content (Ca wt%).Right panel: Relation between and mineral fraction ( m ) as described in Eq. (15)

Fig. 8
Fig. 8 BMDD for an activation frequency Ac.f = 4 BMU/mm 2 /year and mineralisation kinetics T sec equal to A 0.5 years, B 1 year, C 2 years, D 4 years, E 6 years, F 8 years, G 10 years, and H 12 years

Fig. 9
Fig. 9 Comparison of BMU formation and resorption spaces from (left) initial geometry to (middle) R On = 100μ m as in (Martin et al. 1998) and Ac.f = 6 BMU/mm 2 /year and (right) R On = 40μ m and Ac.f = 18 BMU/mm 2 /year as proposed by (Parfitt 1983).Both cases consider T sec = 8 years

Fig. 10
Fig. 10 Bone remodelling simulation using the parameters proposed by Parfitt (Parfitt 1983): R On = 40μ m and Ac.f = 18 BMU/mm 2 /year.A Ca MEAN evolution for different mineralisation kinetics.BMDDs for B T sec = 6 years and C T sec = 8 years

Fig. 12
Fig. 12 Bone microstructure observed for different simulations of osteoporosis: A High FP scenario, corresponding to FP = 475 days (Villanueva et al. 1966).B High Ac.f scenario, corresponding to Ac.f = 8 BMU/mm 2 /year.C Underfilling scenario, corresponding to

Fig. 15
Fig. 15 Short time scales highlight the first mineralisation kinetics.The evolution of the ash factor is little affected by the secondary mineralisation time T sec

Table 1
Discrete bone remodelling model parameters b Local time counter: follows evolution of active BMU m Local time counter: follows the mineralisation phase 1 Local time point: resorption ends / reversal period starts 2 Local time point: reversal period ends / formation startsTable 1 (continued)

Table 2
Input parameters for the discrete bone remodelling model.f.p.: free parameters

Table 3
Mineralisation law parameters for different secondary mineralisation times T sec T sec [years] c 1 [Ca wt%] c max [Ca wt%] t 1 [days] t 2 [days]