Following Pivonka et al. (2008), the BR process is described by a set of differential equations of the cell populations. The BCPM considers the RANK–RANKL–OPG signalling pathway, the action of TGF-β, the mechanobiological feedback on bone cells, the effect of bone mineralisation and the accumulation of microstructural damage by fatigue and its repair through BR. The bone cell types considered in the current model are: osteocytes (Ot), osteoblast precursor cells (Ob_p), active osteoblasts (Ob_a), osteoclast precursor cells (Oc_p) and active osteoclasts (Oc_a). The cell pools of uncommitted osteoblasts (Ob_u) and osteoclasts (Oc_u) are assumed constant. d O b p d t = D O b u ⋅ O b u ⋅ Π act TGF − β + P O b p ⋅ O b p ⋅ Π act ψ b m ⋅ K O b p − D O b p ⋅ O b p ⋅ Π rep TGF − β ⋅ F O b p (1) d O b a d t = D O b p ⋅ O b p ⋅ Π rep TGF − β ⋅ F O b p − A O b a ⋅ O b a (2) d O c p d t = D O c u ⋅ O c u ⋅ Π act RANKL − D O c p ⋅ O c p ⋅ Π act RANKL ⋅ F O c p (3) d O c a d t = D O c p ⋅ O c p ⋅ Π act RANKL ⋅ F O c p − A O c a ⋅ O c a ⋅ Π act TGF − β (4) d O t d t = η d f b m d t (5) where D O b u , D O b p , D O c u and D O c p are the differentiation rates of Ob_u, Ob_p, Oc_u and Oc_p, respectively; while A O b a and A O c a are the apoptosis rate of Ob_a and Oc_a, respectively. η is the concentration of osteocytes in bone matrix, which is assumed constant as in (Martin et al., 2019), thus leading to proportional variations of osteocytes population and fraction of bone matrix volume per total volume, f _bm (Eq. 5). The second term in the right-hand side of Eq. 1 corresponds to proliferation of osteoblast precursors and P O b p gives the maximum proliferation rate. The factor K O b p is introduced following (Buenzli et al., 2012b) and is a saturation function that prevents the proliferation of osteoblast precursors if their population exceeds a certain value O b p max . K O b p = 1 − O b p O b p max i f O b p < O b p max 0 i f O b p ≥ O b p max (6)

The factors Π act TGF − β and Π rep TGF − β represent activator and repressor functions related to the binding of TGF-β to its receptor. Similarly, Π act RANKL is the activator function related to the RANK-RANKL binding, while Π act ψ b m is a function of the mechanical stimulus that regulates the anabolic part of the mechanical feedback in the proliferation term.

The mechanical feedback regulation of bone is included in the model following the Mechanostat Theory proposed by Frost (Frost, 2003) and distinguishing three zones (disuse, no net effect in bone mass and overload) as a function of the “Minimally Effective Strains” (MES). The pathologic overload is indirectly considered through the accumulation of microstructural damage. Mechanical disuse is assumed to enhance the production of RANKL on osteoblast precursors depending on the strain energy density (SED) of bone matrix, through the factor Π act RANKL . Overload is assumed to promote bone formation by proliferation of osteoblast precursors through the activator function Π act ψ b m . The SED, termed here ψ _bm , was used as a measure of the mechanical stimulus sensed by bone cells to drive bone adaptation, as traditionally done in the literature (Huiskes et al., 1987; Beaupré et al., 1990). ψ _bm was used here as an alternative to the strains MES, used in the Mechanostat Theory (Frost, 2003). In a uniaxial stress state both variables are related through: ψ b m = 1 2 E ⋅ M E S 2 (7) E being the Young’s modulus. In a general stress-strain state SEDs is given by: ψ b m = 1 2 σ i j ε i j (8) where the Einstein summation convention has been used and σ _ij and ɛ _ij are the components of the stress and strain tensors, respectively. For brevity, only some of the equations of the previously developed model are shown here. The rest are given in Supplementary Section S3.

The main novelty of the present model is the introduction of the concept of cell availability, a threshold based phenomenon whereby we assume that the differentiation of precursors into mature cells is activated only when the population of the former reaches a certain upper threshold and continues until that population falls below a lower threshold. Thus, two binary variables F O b p and F O c p are defined to activate (F_X = 1, with X = Oc_p or Ob_p) or deactivate (F_X = 0) the differentiation of osteoblasts and osteoclast precursor cells, respectively. if O b p > O b p u p , t h → F O b p = 1 until O b p < O b p low,th → F O b p = 0 (9a) if O c p > O c p u p , t h → F O c p = 1 until O c p < O c p low,th → F O c p = 0 (9b)

We note that in Eq. 6 O b p max > O b p u p , t h . Therefore, following Eqs 9a, 9b, precursor cells accumulate until their population reaches the upper threshold, when they begin to differentiate into their respective mature (and active) forms. Then, the precursor cell pool drops to the lower threshold, when the concentration is not sufficient to ensure the continuity of the process and at that moment differentiation stops. At this point, the concentration of precursor cells can rise again if the necessary environmental factors are met and a new remodelling cycle can start if the upper threshold is reached again.

The action of mature osteoclasts and osteoblasts in the BMU is reflected by temporary or permanent changes in porosity. Bone matrix fraction is defined as the volume of bone matrix, V _b , per total volume of the bone sample, V _T , expressed as a percentage, i.e., f b m ( % ) = V b V T ⋅ 100 . Its variation is obtained through the balance between resorbed and formed tissue: d f b m d t = − k res ⋅ O c a + k form ⋅ O b a (10) where k _res and k _form are, respectively, the rates of bone resorption and osteoid formation. Figure 1 shows a schematic representation of the BCPM model (Martínez-Reina et al., 2021b) on which the current spatio-temporal model is based. The values of the constants of the BCPM are given in Supplementary Table S1.

The mathematical framework described above can be applied in a bone RVE to reproduce the different phases of an isolated BMU acting locally within the RVE. However, if the RVE is large enough to hold several BMUs, modelling the discrete activation of a single BMU is no longer meaningful and, more importantly, the interactions between those BMUs cannot be modelled in the sense addressed in this work. For this reason, a 3D domain must be considered to account for these interactions. In particular, a 3D FE mesh is proposed in a non-local approach, with each element being small enough to hold at most one BMU. In this way, the movement of the BMU can also be modelled by changing the element that contains the BMU at each time.

In a FE mesh, the conditions for activation, deactivation and progression of a BMU should be determined. In the local approach, addressed in the previous section, the only condition to be fulfilled in order to activate a BMU was that the concentration of Oc_p must be higher than a certain threshold for them to differentiate into Oc_a. In the following, an element [or better an integration point (IP)] ¹ where the differentiation of precursor cells into mature cells is taking place will be called “active” and this includes active resorption and active formation, since each process is evaluated separately. For the non-local approach proposed in this work, other conditions are imposed in order to prevent the activation of a certain process (resorption or formation) if the same process is active in a nearby element. We hypothesize that a new BMU will not be activated if an existing BMU is progressing in the neighbourhood. Instead, the metabolic cost of activating a new BMU will be employed in “feeding” the existing BMU, with the recruitment of more progenitor cells to the resorbing and formation sites, respectively. To do so, the model must take into account what is occurring in the vicinity of a given IP and not only the local events.

The neighbourhood is defined following the algorithm proposed by Calvo-Gallego et al. (2021). For a given IP, p, a neighbourhood of IPs, NBH ^p , is defined by the IPs contained in a sphere of radius R centered at p: N B H p = i ∈ IPs of FE mesh , i = 1 , … , N p ; such that ‖ r i − r p ‖ ≤ R (11) where r _i denotes the position of IP i. Next, the resorption and formation fronts of the BMU are distinguished. To this end, two state variables are defined at each IP: resorption state, RS, and formation state, FS, so that: R S = 1 if this IP is the resorption front ⇒ F O c p = 1 2 resorption is active but this IP is not the front ⇒ F O c p = 1 0 resorption is inactive at this IP ⇒ F O c p = 0 (12) F S = 1 if this IP is the formation front ⇒ F O b p = 1 2 formation is active but this IP is not the front ⇒ F O b p = 1 0 formation is inactive at this IP ⇒ F O b p = 0 (13)

The procedure to establish if one process (either resorption or formation) is activated or deactivated at a certain IP at the instant t _a is described next. It is important to note that resorption and formation are treated separately. For this reason, in the following procedure, the state variables are designated as XS, where X can stand either for R or F. The coupling of both processes and the correct sequence (resorption, reversal, formation) is not enforced, but emerges as a model output.

1. First, the IPs of the mesh are ordered from the highest to the lowest concentration of precursors Oy_p, where y stands for b or c depending on the process being activated. In the next steps the IPs will be analysed in this order. The rationale for this is explained in detail in Figure 2.

Next, the following conditions are checked for the subset C.

3. The concentration of precursor cells Oy_p at IP j is the highest in NBH ^j .

These are the conditions for the activation of a new process. Once activated, the progression of each process along the mesh is subjected to the fulfillment of the following criteria, otherwise the front will not move.

5. The concentration of precursor cells (osteoclasts or osteoblasts) of the destination IP must be above the corresponding upper differentiation threshold.

If these criteria are fulfilled, the destination IPs become the new fronts (XS = 1) and the origin IP becomes a simple active point XS = 2. It must be noted that a smooth progression of the BMU is not guaranteed. Indeed, the BMU could theoretically leap one or more IPs if the radius of the neighbourhood spans more than 2 elements and this would be a limitation of the model which can be avoided by choosing a sufficiently small radius R.

9. A certain process is deactivated at an IP if the concentration of precursor cells (osteoclasts or osteoblasts) falls below the lower differentiation threshold. In this case XS changes from 1 or 2 to 0.

The activation of a BMU starts when the first osteoclasts are differentiated from their precursors (i.e., when the upper threshold O c p u p , t h is exceeded and F O c p = 1 ) at an IP of a neighbourhood where there was no other existing BMU. At this IP RS becomes equal to 1.

To summarize the algorithm, each process of the BMU (resorption or formation) is activated at the point where the concentration of precursors is the highest to proceed with the differentiation of precursors into active cells. This differentiation will proceed along the gradient of precursors concentration until it is deactivated, when the concentration of precursors is below the lower threshold. In view of Eq. 3, an increase of RANKL production will lead to a rise in the concentration of osteoclast precursors. One possible reason for that RANKL increase is the accumulation of microstructural damage (Martínez-Reina et al., 2021b) through the factor Π act dam (see Supplementary Section S2.3 for more details). Thus, one of the reasons for the BMU resorption front to move in this model is to repair highly damaged areas, in accordance with the theory of targeted bone remodelling (Parfitt, 2002). On the other hand, the concentration of osteoblast precursors will increase in overloaded areas, as the proliferation term in Eq. 1 is upregulated by mechanotransduction through the factor Π act ψ b m . Ob_p will also increase in areas where resorption has recently taken place, as osteoclasts release TGF-β from bone matrix and this upregulates the differentiation of uncommitted osteoblasts into osteoblast precursor cells through Π act TGF−β (see Eq. 1). Hence, the BMU formation front will follow the path previously set by osteoclasts until formation is deactivated, so allowing the different phases in sequential order: activation, resorption, reversal, formation and quiescence.

The role of TGF-β in this sequence is paramount as will be shown later in the results. Its concentration in the bone compartment is governed by the following differential equation in the original model (Pivonka et al., 2008): d T G F − β d t = P TGF − β + α TGF − β ⋅ k res ⋅ O c a + D ̃ TGF − β ⋅ T G F − β (14) where P ^TGF−β is the external production of TGF-β. The second term on the right-hand side represents the release of TGF-β from bone matrix through bone resorption, with α _TGF−β being the concentration of TGF−β in bone matrix. Finally, D ̃ TGF − β is the TGF-β degradation rate. In temporal-only BCPMs, stationarity of Eq. 14 is assumed in order to compute the TGF-β concentration in the RVE. ³ However, assuming stationarity in the current model would neglect the effect that the release, action and degradation of TGF-β has on the coupling of the BMU processes. Therefore, the differential Eq. 14 is solved in the current model. We note that Buenzli et al. (2012a) also pointed out that the distribution of TGF-β from the BMU cutting cone to the closing cone has a significant effect on bone cell density distribution and separation of cell populations.

The new spatio-temporal BCPM model was tested in a 3D domain (see Figure 3). It was developed in ABAQUS and consists of a parallelepiped of dimensions 0.9 × 0.9 × 1.8 mm³ (see Figure 3). It was meshed with 6,750 eight-noded hexahedral isoparametric elements, with reduced integration (one integration point per element, named C3D8R in the ABAQUS Element Library). All hexahedra are regular of side 60 μm. The choice of this size was based on the fact that a BMU progresses at a rate between 20 μm/day (Parfitt, 2002) and 40 μm/day (Martin et al., 2015). Besides, the model corresponds to a piece of trabecular bone which is analysed at the mesoscale, i.e., from a Continuum Mechanics point of view and without considering its microstructure, though accounting for the spatial variation of f _bm across the trabecular structure at the mesoscale. For that purpose, an element size of 60 μm is also adequate.

Three load cases were analysed: uniaxial compression ( σ z z 1 = − 0.5 MPa), uniaxial tension ( σ z z 1 = + 0.5 MPa) and compression plus bending ( σ z z 1 = − 0.5 MPa and σ z z 2 = − 1.33 N/mm, so that the longitudinal stresses ranged from −0.2 to −0.8 MPa). The boundary conditions were chosen to ensure a uniform uniaxial stress state in the first two cases. More precisely, the Z displacement of all nodes of one face (red dots, see Figure 3 left); the Z and Y displacements in one corner (blue dot) and the Z and X displacements in the opposite corner (black dot) were restrained to prevent the rigid-body motions. The same boundary conditions were applied to the compression plus bending case.

A subregion of this model (blue in Figure 3 right), far enough from boundary conditions and concentrated loads, was selected to evaluate the results. This subregion is a parallelepiped composed of 504 elements (6 × 6 × 14) and placed in the core of the model. This subregion is hereinafter referred to as the domain under study. Two cases have been simulated: a homeostasis situation and a model including a highly damaged region.

The thresholds in Eq. 9 were calibrated to obtain a BMU activation frequency and lengths of the BMU phases in accordance with the literature. For this purpose, the thresholds were varied in the following ranges: O b p u p , t h ∈ 0.04 , 0.02 , O b p low,th ∈ 0.0025 , 0.0005 , O c p u p , t h ∈ 0.025 , 0.005 and O c p low,th ∈ 0.0025 , 0.0005 . These ranges were divided in 3 equally spaced subranges and all the combinations of the resulting 4 points of each threshold were evaluated. Subsequently, a local gradient search was performed starting from the best combination, being the total difference with the lengths of Table 2 the objective function to be minimised in this search. The optimum is shown in Table 1 and this was used for the nominal values of the thresholds. A sensitivity analysis of the thresholds was performed around those nominal values to evaluate their effect on the length of the phases.

To check whether homeostasis was reached after the initial perturbation, the evolution of f _bm was analysed in the domain under study for the compression load case. First, the average of f _bm was calculated for the elements of the domain under study to yield f b m aver ( t ) . The temporal evolution of f b m aver ( t ) showed an initial adaptation to the applied load, which was low for the initial distribution of f _bm , thus resulting in a decrease of bone mass. During this phase of adaptation very pronounced fluctuations were observed as a consequence of the numerous remodelling events occurring within the domain. Thus, a moving average filter was applied to the temporal evolution of f b m aver ( t ) to give: f b m m v t = ∑ t − n / 2 t + n / 2 f b m aver t (15) where n is the number of time points considered for the average (or time window). Note that t + n/2 cannot exceed the simulation time (t _max ) and t − n/2 cannot be less than 0 and therefore, the window must be truncated at those endpoints.

The evolution of ψ _bm was also obtained in the domain under study, to analyse the relationship between f _bm and the mechanical stimulus. ψ b m aver and ψ b m m v were calculated through spatial and temporal averaging of ψ _bm , analogously to f b m aver and f b m m v . Figure 5A shows the temporal evolution of f b m m v after the filter was applied for two different time windows. Figure 5B compares f b m m v ( t ) and ψ b m m v ( t ) . Figure 5 corresponds to the compression load case with an initial random distribution of f _bm ∈ [10, 80]%. The evolution of f b m m v ( t ) when initial f _bm ∈ [27, 29]% is shown in Figure 6 for comparison.

Figure 7 represents the temporal evolution of resorbed bone volume and formed bone volume per unit volume and unit time (i.e., RBV = k _res ⋅Oc_a and FBV = k _form ⋅Ob_a), for an element inside the domain under study.

Several activations of BMUs took place in the domain under study, though not simultaneously. The average length of the resorption, inversion (or reversal), formation and quiescence phases was calculated from the results of all the elements inside the domain under study. The first phase corresponds to a transient period during which the bone volume fraction is adapting to the applied load and must be ruled out to focus on the homeostatic situation, which was assumed to be reached around day 10,000 (see Figure 5). The average lengths of the phases were calculated from day 13,000 until the end of the simulation (day 15,000). These averages are given in Table 2 and compared with experimental data from the literature. The activation frequency is calculated in histomorphometric studies from the average length of the BMU phases (Eriksen et al., 1990; Steiniche et al., 1994) and represents the number of BMUs passing through a given bone site per year: f act,hist B M U s y e a r = 365 T R + T I + T F + T Q (16) where T _R , T _I , T _F and T _Q are the average resorption, inversion (reversal), formation and quiescence periods. The values obtained for this frequency are given and compared with the literature in Table 2.

The 3D BMU activation frequency, f _act,3D , was calculated by counting the number of activations of BMUs occurring within the domain under study per unit time and unit volume. This 3D activation frequency can be converted into a 2D activation frequency following Martin (1984) who related both through: f act,3D = k ⋅ f act,2D (17) with k being a length parameter. Later, Hernandez et al. (1999) used k = 1, so identifying the 2D and 3D frequencies. If this assumption is accepted, the estimated values of f _act,2D would be in the range given by Frost (1964) (see Table 2).

The results of the sensitivity analysis are presented in Table 3 for the case of compression. The influence of the thresholds and R (the radius of the sphere defining the neighbourhood) was studied.

Regarding the simulation of the case of absence of TGF-β, the evolution of the resorbed bone volume per unit time and unit volume (RBV) and the corresponding formed bone volume (FBV) is plotted in Figure 8 for an element inside the domain under study. Compared to Figure 7, where the cycles had a reasonable shape and followed the normal sequence resorption—inversion—formation—quiescence, now the BMU cycles appeared totally uncoupled. Different rare events can be seen in Figure 8: several formation cycles taking place (before or after a resorption cycle), resorption and formation occurring simultaneously and cycles with highly variable resorption and formation lengths.

To quantify the occurrence of these rare events, three types of resorption/formation cycles have been defined.

• A BMU cycle (or remodelling event) is regarded normal if a resorption cycle is followed by a formation cycle in less than T I max = 20 days, the maximum admissible inversion time.

All these cycles have been counted in each element of the domain under study in the normal case and in the case of absence of TGF-β. The results are presented in Table 4 for the three load cases. The temporal evolution of f _bm in the absence of TGF-β (not shown) is similar to that of Figure 5A, though with more irregular oscillations.

In order to analyse how the BMU progresses along bone matrix repairing the damaged tissue, we have studied the BMU front and those elements where the differentiation from precursor cells to mature cells is active. In Figure 9, the line of highly damaged elements are shown at different times during the resorption cycle: in blue, those where differentiation is not active; in red, the BMU front, with active differentiation; in green, those elements with active differentiation but which are no longer the BMU front. It can be noted that the BMU was activated at the element with the highest damage (recall Figure 4). Later, the BMU front progressed along the damage gradient, towards the next elements with high damage. Differentiation of mature osteoclasts keeps active in those elements through which the BMU front runs, so resorbing bone and repairing damage until the precursors concentration falls below the lower threshold, when they become inactive.

The evolution of the damage variable over time in the compression load case can be seen for the damaged elements in Figure 10. It can be noted how a fraction of damage is repaired in the first element when the BMU is activated and how the repair process progresses along the line of damaged elements. Between 20% and 25% of the original damage is repaired.

Our results demonstrated that the spatio-temporal BCPM model is able to reach a homeostatic state in a 3D spatial domain, characterised by a slightly oscillating variation of f _bm which remains constant on average (see Figure 5A) due to an equilibrium between the resorption and the formation cycles of BMUs. This is opposed to the monotonic and smooth way the previous temporal-only model reaches a homeostatic state after a parameter perturbation. The oscillations are more realistic since bone turnover events occur sporadically and spread out across bone sites, as was demonstrated in our model simulations, but the amplitude of the oscillations diminishes with the size of the domain where the results are averaged.

It is important to note that during the first phase of the simulation shown in Figure 5, the average f _bm is adapting to the external load, which is too low for the initial distribution assumed for f _bm . This leads to a global bone loss during a relatively long transient period where fluctuations of f _bm are pronounced. These fluctuations are damped out in the long-term, showing that the model is able to reach an equilibrium state adapted to the applied load. The amplitude of the oscillations and the length of the transient period depend on the applied load compared to the homeostatic load and on the initial distribution of the variables. In order to consider a situation more realistic than the uniform distribution, some variables were randomly distributed in space. In the case of damage and concentration of osteoclast precursors the amplitude of this perturbation was based on the normal values achieved with the previous BCPM model (Martínez-Reina et al., 2021b). The perturbation of damage was quickly damped, having no effect in the mid-term. The perturbation of other variables such as mineral content also damp out very quickly, so they were not included in these simulations. The perturbations of Oc_p and specially of f _bm have a stronger effect since they make BMU activations be more uniformly distributed over time throughout the domain under study. The transient period will be left out from the rest of the discussion and we will focus on the equilibrium state.

Starting from an almost uniform distribution of f _bm (in the range [27, 29]% in Figure 6) leads to an almost synchronous activation of BMUs with large fluctuations of f _bm , which is not realistic. This initial distribution of f _bm is also unrealistic from another point of view, as it cannot correctly represent the spatial variation of porosity throughout the trabecular microstructure with elements as small as those of 60 μm in size. In conclusion, an appropriate initial distribution of f _bm is crucial to obtain realistic results and should be chosen according to the element size.

The proposed algorithm sorts the IPs based on the precursors concentration, Oy_p. As explained in Figure 2 a gradient of Oc_p could prevent the activation of BMUs in a spurious way. For instance, consider a monotonic gradient of Oc_p in z-direction of the model (see Figure 3). If the order of calculation of IPs started in the growing direction of Oc_p, BMUs would only be activated at the end section and therefore, the activation frequency would be highly underestimated, though in a spurious way. This monotonic gradient of Oc_p is not realistic and in fact it is difficult to achieve with the random assignment of initial Oc_p. With this random assignment, the activation frequency is just slightly underestimated, but this situation must be prevented in any case and for that reason sorting the IPs in the algorithm is important. Only in this case, BMU activation is allowed or prevented depending on the concentration of precursors, in the IP and the surroundings, and not influenced by the order of calculation of IPs.

Figure 5B allowed to compare the temporal evolution of f _bm and the stimulus ψ _bm . It could be seen how the minima of f _bm coincide with the maxima of ψ _bm . As the resorption phase progresses and f _bm decreases in a certain region, the remaining tissue becomes overloaded and the mechanical stimulus (SED) rises, thus promoting the proliferation of osteoblast precursors and eventually the onset of the formation phase. During this phase that trend is reversed as the rise of f _bm increases the stiffness, so reducing the strains and consequently the SED (note that the stress is kept constant). During the quiescence phase (f _bm ∼ constant) the slow fall of ψ _bm follows the gradual increase in stiffness caused by the mineralisation of the newly formed tissue.

The different phases of a BMU are simulated with the current model in a realistic manner. First, the BMU is activated at a certain location with the differentiation of osteoclast precursors into mature osteoclasts, so starting a resorption cycle aimed at repairing the damaged tissue. The highest peaks observed in Figure 7 correspond to the resorption cycles, which are higher and much shorter than the formation cycles, so revealing that resorption is faster than formation but it extends over a shorter period of time. The new model proved also valid to simulate the reversal phase, so separating the resorption and formation events. Thus, once osteoclasts have finished the resorption phase at a site, they move or undergo apoptosis and the reversal phase takes place before the formation phase begins. Finally, the formation cycle ends and the remodelled region remains in a quiescence state for quite a long time (around 1.5–2 years) until a new BMU is activated. The areas under the RBV and FBV curves are equal, so indicating that the resorbed volume in a cycle is equal to the formed volume (homeostasis). Moreover, the thresholds introduced in the BCPM model were sufficient to obtain a set of phase lengths similar to those found in the literature (see Table 2).

As stated in the introduction, the reversal phase can last from 8–9 days (Eriksen et al., 1984; Hazelwood et al., 2001) to several weeks (Delaisse, 2014; Martin, 2021). The discrepancy in these experimental results may be due to a different interpretation of the reversal phase, with the latter studies including in it the mineralisation lag time. Similar discrepancies exist also in the length of the formation period, with values ranging from 64 (Hazelwood et al., 2001) to 174 days (Eriksen et al., 1990). Even in the latter study the range between the 10% and the 90% percentiles spans 74–481 days. Motivated by these discrepancies, we have used here the values provided by Hazelwood et al. (2001) for comparison of the average times obtained for the resorption, formation and reversal periods, as these authors compared many histomorphometric studies.

Some histomorphometric studies calculated the BMU activation frequency as the inverse of the total time of the BMU cycle. i.e., the summation of the resorption, reversal, formation and quiescence periods (Eriksen et al., 1990; Steiniche et al., 1994). Thus, as the results obtained in this work for those periods are similar to those given the literature, so is the activation frequency measured in this manner. The 3D activation frequency was also calculated here, although it is not possible to compare it with any experimental results. Taking into account the relationship between the 2D and the 3D activation frequencies proposed by Hernandez et al. (1999), our result would be in agreement with the values provided by Frost (1964). This author measured the evolution of BMUs activation rate with age, obtaining values between 1 and 2 BMUs/year/mm² for people over 30 years old. The mean value obtained in this work is in that range; however, it is important to note that the relationship between the 2D and 3D frequencies is not absolutely reliable because the parameters were estimated with a high degree of uncertainty, as stated by the authors (Martin, 1984; Hernandez et al., 1999).

The sensitivity analysis whose results were given in Table 3 shows that O c p u p , t h is the threshold with the greatest sensitivity. It affects significantly the inversion and quiescence times and to a lesser extent the formation time and the activation frequency. It can be noted that O b p low,th has only a slight influence on the average phase lengths, as does O c p low,th , except maybe for the reversal time. O b p u p , t h has a high influence on the reversal time, while the rest of times and the activation frequency remain roughly unaffected, except if O b p u p , t h is increased by 20% or more. It could be checked that the abnormally high values of T _I were caused by the same lack of coordination of the BMU cycles observed in the simulations with no TGF-β. The radius of the sphere defining the neighbourhood, R, has a strong influence on the 3D activation frequency: the larger the radius, the smaller the 3D activation frequency. Certainly, as R increases, new activations are prevented in a larger volume if one BMU is already active within the neighbourhood. Finally, both the thresholds and R seem to have a small influence on the length of the resorption period.

The influence of R on the activation frequency is related to the features of the FE mesh, particularly to the element size and the integration scheme (reduced vs. full integration). If these are changed, the suitable R may also change. The reason for this is that both variables also affect the IPs contained within the neighbourhood. In fact, the most influential factor on the activation frequency is the distance from the centre of the neighbourhood to the first IP outside it, where the BMU activation is no longer prevented.

TGF-β is a cytokine stored in the bone matrix and released by osteoclasts during bone resorption, which has been reported to play a crucial role in the coupling formation to resorption (Durdan et al., 2022). It is known that this cytokine promotes osteoclast apoptosis and differentiation of uncommitted osteoblast into osteoblast precursors, while it downregulates the differentiation of osteoblast precursors into mature osteoblasts (Pivonka et al., 2008). In other words, TGF-β promotes the accumulation of osteoblast precursors, so preparing the upcoming formation cycle. The upregulation of osteoclasts apoptosis also helps bring the ongoing resorption cycle to an end. Therefore, TGF-β seems to be of paramount importance in the synchronization of the BMU cycle. We have shown with our model that in the presence of TGF-β BMU cycles exhibit a normal pattern, with resorption followed by formation after a reversal phase and a long quiescence period before the next resorption cycle (see Figure 7). Moreover, the average time for each phase is in agreement with the literature. However, in the absence of TGF-β (see Figure 8), the BMU cycles become completely uncoupled. Resorption and formation would take place simultaneously at a given bone site, which would not have any biological sense. The reversal phase could be extremely long and even formation or resorption cycles could appear isolated. Moreover, these rare events are not one-off. On the contrary, they occur very frequently, the normal cycles being the less usual ones, as can be seen in Table 4.

The comparison of Figures 7, 8 might lead one to think that the equilibrium of f _bm is not achieved in the absence of TGF-β, as the areas under the FBV and RBV curves are not equal. However, it should be noted that the evolution of Figure 7 is not representative of what occurs in average, but only of this particular element and time window. If f _bm is averaged in the domain under study over time, an approximately constant f b m m v is also obtained in the equilibrium, though with larger and more chaotic oscillations than those observed in Figure 5.

As stated before, a lack of coordination of the BMU cycles was also observed in some cases of the sensitivity analysis where T _I was abnormally high (see Table 3). However, the percentage of normal cycles was much higher in these cases and more importantly, that lack of coordination can be avoided with a proper choice of the thresholds, something that is not possible when TGF-β is absent.

Another phenomenon that is better described in intermittent models is the absorption of bisphosphonates. These drugs are also stored in the bone matrix and released through bone resportion, like TGF-β, though their effect is different and seems to affect only osteoclast by impairing its resorptive capacity (Halasy-Nagy et al., 2001; Takagi et al., 2021) and enhancing its apoptosis rate (Sato et al., 1991; Halasy-Nagy et al., 2001; Yu et al., 2018; Takagi et al., 2021). In a recent work we have developed a pharmacokinetics-pharmacodynamics (PK-PD) model to simulate the absorption of alendronate and its effect on bone turnover. This model implements a continuous BCPM and therefore, it does not allow to simulate well-known features of alendronate such as its higher deposition rate at sites where bone turnover is more intense (Porras et al., 1999). This is feasible in spatio-temporal models such the one proposed here as remodelling sites are easily distinguished from quiescent sites. This would lead to a heterogeneous distribution of the drug that could also allow to consider other important effects, such as the saturation of the drug absorption in certain sites, which is not feasible in temporal-only BCPM. Analogously, the heterogeneous distribution of mineral within bone matrix also emerges from the specific localisation of the remodelling events and therefore, it can be addressed by spatio-temporal bone remodelling models, but hardly by temporal-only models.

We have simulated the activation, progression and deactivation of BMUs in a 3D FE model of a piece of trabecular bone with f _bm ∼ 30% under different load cases. This bone volume fraction was selected based on the experimental data taken to compare the quiescence time that was measured in cancellous bone (Steiniche et al., 1994). Trabecular thickness is approximately 90 μm (Vesterby et al., 1991) and the element size of the FE mesh is 60 μm. Therefore, modelling the porous bone as a continuum allows the BMUs to move in any direction, whereas in reality they could only progress onto the existing bone tissue, being forced to turn to follow the trabecula, something that does not occur in our model, which does not account for the microstructure of the trabecular bone. This requires a more complex model with geometrical constraints to the movement of the BMUs and is left for future works.

The BCPM presented in this work is able to model the activation of BMUs in highly damaged areas and to steer them to follow a crack (or a region with high damage) in order to repair that damage, as the theory of targeted bone remodelling holds. It is important to note that not several but only one BMU was activated to follow the crack path. The model was designed to avoid the multiple activation of BMUs in a small region, which would not be meaningful from a biological or metabolic point of view. This means that osteoclasts continue to be recruited to the existing BMU if the mechanobiological stimulus (the presence of damage in this case) is still high enough.

As presented in the results of the crack simulation, the BMU repairs only 20%–25% of the existing damage. This result is not very realistic from a biomechanical point of view and this constitutes a limitation of the study. Indeed, osteoclasts dissolve bone mineral by secreting acids and degrade the organic matrix with specialized proteinases (Blair, 1998). If the resorption front spanned the entire damaged region, the resulting resorption cavity would be a void within the mesh the size of several elements, similar to the voids that may appear implementing the element killing technique. This can lead to numerical issues in our FE simulations, arising from the zero stiffness implied by the total removal of bone tissue and was disregarded here. Our model cannot contemplate resorption cavities as such and treats porosity in the Continuum Mechanics sense, as a variable that represents the volume occupied by pores within the elements, without representing the pores explicitly. Similar to that, damage is also considered in the Continuum Mechanics sense, as a measure of microcracks density. Such model is unable to repair damage completely and for this reason, it would be more appropriate to treat this case as the remodelling of a region of high diffuse damage rather than a crack.