A bone remodelling model including the effect of damage on the steering of BMUs

A bone remodelling model including the effect of damage on the steering of BMUs J.Mart´ınez-Reinaa,∗, I. Reinaa, J. Dom´ıngueza, J.M. Garc´ıa-Aznarb aDepartment of Mechanical Engineering, Universidad de Sevilla. Escuela T´ecnica Superior de Ingenier´ıa, Camino de los Descubrimientos s/n E-41092 Sevilla, Spain bDepartment of Mechanical Engineering, University of Zaragoza, Arag´on Institute of Engineering Research (I3A) Abstract Bone remodelling in cortical bone is performed by the so-called basic multicellular units (BMUs), which produce osteons after completing the remodelling sequence. Burger et al. (2003) hypothesized that BMUs follow the direction of the prevalent local stress in the bone. More recently, Martin (2007) has shown that BMUs must be somehow guided by microstructural damage as well. The interaction of both variables, strain and damage, in the guidance of BMUs has been incorporated into a bone remodelling model for cortical bone. This model accounts for variations in porosity, anisotropy and damage level. The bone remodelling model has been applied to a finite element model of the diaphysis of a human femur. The trajectories of the BMUs have been analysed throughout the diaphysis and compared with the orientation of osteons measured experimentally. Some interesting observations, like the typical fan arrangement of osteons near the periosteum, can be explained with the proposed remodelling model. Moreover, the efficiency of BMUs in damage repairing has been shown to be greater if BMUs are guided by damage. Keywords: Bone remodelling, basic multicellular units, cortical bone, orientation of osteons ∗Corresponding author Email address: [email protected] (J.Mart´ınez-Reina) Preprint submitted to Journal of the Mechanical Behavior of Biomedical MaterialsDecember 27, 2013 1. Introduction. Simulation of bone remodelling is a classic problem in Biomechanics but it continues receiving attention (see e.g. (Hartmann et al., 2011; Kaczmarczyk et al., 2011; Malachanne et al., 2012)) for its importance in clinical applications (see e.g. (Prendergast et al., 2011; Tomaszewski et al., 2012)) and because many questions are still open and need further investigation. In this work some of those questions are addressed to propose a model of bone remodelling in cortical bone. Cortical bone is usually assumed to be a transversely isotropic material, being its mechanical properties strongly influenced by the orientation of osteons. More precisely, the longitudinal axis of the osteons is the perpendicular to the plane of isotropy (Yoon and Cowin, 2008a, 2008b). The orientation of an osteon is given by the trajectory followed by osteoclasts and osteoblasts, which act sequentially in an association of cells called BMU (Basic Multicellular Unit) to form the osteon. This paper provides some insights about some mechanisms that could regulate that trajectory. The direction followed by BMUs in their progression is not clear. Burger et al. (2003) proposed that they follow the direction of the prevalent local stress: either maximum or minimum principal stress, whichever has a higher absolute value. Their argument is based on the fact that osteoclasts are attracted to sites of low mechanical stimulus, with osteocytes playing an important role in this attraction. If the trajectory of the BMU is aligned with the maximum stress direction (or minimum stress direction in compressive situations), the area in front of the cutting cone is a low stimulus area (recall the stress distribution around a hole subjected to far-field stresses in one direction) and, thus, more osteoclasts will be attracted to the front of the cutting cone to continue resorbing bone in the maximum (or minimum) principal stress direction. This alignment of osteons with principal stress directions has been implemented in some bone remodelling models (van Oers et al., 2008a, 2008b) including one developed by the authors (Mart´ınez-Reina et al., 2009). In the case of long bones, compressive, bending and torsional stresses would lead, on that assumption, to trajectories that spiral slightly around and through the cortex at low angles in opposite directions on the lateral and contralateral sides of the bone, merging to a longitudinal trajectory anteriorly and posteriorly (Heˇrt et al., 1994; Petrt´yl et al., 1996). However, Martin (2007) showed that BMUs must be driven by another equally important factor: damage. More precisely, with a clear objective: 2 damage removal. This author claimed that, in case BMUs tunnelled through cortical bone at random locations, but following straight trajectories, they would need to have a greater cross-section (40 times) than they actually have, in order to repair damage with the efficiency observed experimentally. Another plausible explanation for that high efficiency would be, following Martin, that BMUs do not resorb bone at random locations, but they are steered toward the damaged areas. In this guidance, apoptotic osteocytes would play a key role, as they did in the hypothesis of Burger et al. (2003). In this case, osteocytes close to the damaged areas might undergo apoptosis and produce biochemical signals that attract osteoclasts, this way steering the BMU toward the damaged area. In this work, we propose a new bone remodelling model (from this point on referred to as strain and damage-based model) incorporating a possible steering of BMUs by damage (apart from by principal strain directions) into a previously developed bone remodelling model, referred to as old model (Mart´ınez-Reina et al., 2009). Other modifications have also been included in the strain and damage-based model and are explained in detail in the next section. Moreover, a variation of the strain and damage-based model, named strainbased model and obtained by deactivating the steering of BMUs by damage, is also presented. Both, the strain and damage-based model and the strainbased model have been applied to the diaphysis of a human femur subjected to normal walking loads. The general objective of this paper is to analyse how the guidance of BMUs regulated by damage may explain the actual orientations of osteons in long bones. In particular, to show how this guidance can explain certain phenomena that the hypothesis of Burger et al. (2003), of osteoclasts following the direction of the prevalent local stress, cannot. These phenomena are: the high efficiency of BMUs in damage repairing, confirmed by the observations of Martin (2007), and the typical fan arrangement of osteons near the periosteum. 2. Materials and Methods. 2.1. Old model. Activation and remodelling effect of BMUs The old model (Mart´ınez-Reina et al., 2009) is an extension of a previous isotropic one (Garc´ıa-Aznar et al., 2005) to the anisotropic case. Such extension was intended to include the orientation of osteons and its influence on the 3 anisotropy of bone tissue. The activation and the remodelling effect of BMUs is similarly implemented in both models (isotropic and anisotropic) and is briefly explained next. Nevertheless, consulting both papers (Mart´ınez-Reina et al., 2009; Garc´ıa-Aznar et al., 2005) is advised for a further comprehension of the concepts presented in this section. Bone remodelling is carried out by the so-called BMU (Basic Multicellular Unit), a temporal association of tissue-resorbing osteoclasts and tissueforming osteoblasts working in a coupled and sequential way. The number of active BMUs is controlled by their lifespan, σL, and the BMU activation frequency, ˙ NBMU . This frequency provides the number of BMUs activated per unit time and unit volume and is given in the old model by: ˙ NBMU =fbio(1 −s)Sv(1) where Svrepresents the bone specific surface (bone matrix surface per unit volume). This variable is a function of porosity (Martin, 1984) and its inclusion in equation (1) accounts for the fact that the activation of BMUs takes place on the bone matrix surface. fbio is a parameter including all the biological and metabolic factors that influence on the activation of BMUs and is assumed to be a constant in this work. Finally, srepresents the intensity of a certain mechanobiological signal which inhibits the activation of BMUs. This model follows the inhibitory theory proposed by Martin (Martin, 2000), which establishes that bone lining cells are inclined to activate BMUs and therefore initiate bone remodelling, except in case they receive an inhibitory signal emitted by osteocytes and transmitted through the canalicular network. This signal is defined here as s=ξ ξ+c(1 −diso)a(2) and depends on the mechanical stimulus, ξ∈[0,+∞), and the level of microstructural damage, diso ∈[0,1], which is a measure of the density of microcracks and is related to the loss of stiffness (see (Mart´ınez-Reina et al., 2009) and (Garc´ıa-Aznar et al., 2005)). The constant cis identified with a reference mechanical stimulus, ξ∗, and is positive constant, like the constant a. Thus, s∈[0,1]. In normal conditions, osteocytes transmit the inhibitory signal to the lining cells and BMUs are not activated. However, if they undergo apoptosis, either by damage (dhigh) or disuse (ξlow), the inhibitory signal is low, resulting in a high activation frequency. So, a sufficient number 4 of BMUs will be activated in order to repair the microstructural damage or to resorb the excess of bone, respectively. That inhibitory mechanism is not very clear and, therefore, it makes more sense to work with an activating signal rather than with an inhibitory signal. That activating signal,s∈[0,1], has been defined as: s= 1 −s(3) Following Mikic and Carter (1995), the stimulus ξis defined as the daily strain history, based on the strain level, ¯εi, and the number of cycles, ni, of each load case i: ξ=(∑ i ni¯εm i)1/m (4) The parameter mis taken to be 4 (Whalen et al., 1988) and the effective strain, ¯εi, is defined as a function of the strain energy density, Ui, and the Young’s modulus, E, as follows: ¯εi=√2Ui E(5) The cells involved in the activity of BMUs proceed in a strict order at a certain bone location: after the osteoclasts are recruited, they begin to resorb the old bone tissue, progressing in two directions: radial and longitudinal (if the osteon is approximated to a cylinder). At a given cross-section of the osteon the resorption phase spans a time interval of length TR. Then, after a reversal time of length TI, the formation phase begins, spanning a time TF. This sequence was originally termed as A-R-F sequence (activation, resorption, formation) by Frost (1964a) and lasts TR+TI+TFuntil the crosssection is completely remodelled. While osteoblasts are forming new tissue and closing the resorption cavity at a certain section, osteoclasts continue progressing, until their apoptosis. They progress along a given direction at a rate VBMU during σLdays and always followed by the osteoblasts, which so refill the resorption cavity. The volume of tissue resorbed and formed at time tper unit time and unit volume is given, respectively, by equations (6a) and (6b): ˙vR(t) = ∫t 0 (1 −p)˙ NBMU (t′)AR(t”) VBMU fcdt′(6a) ˙vF(t) = ∫t 0 ˙ NBMU (t′)AF(t”) VBMU fbdt′(6b) 5 where t” = t−t′is the time elapsed since the activation of the BMU. fcand fbare two parameters measuring the activity of osteoclasts and osteoblasts, respectively. It is assumed that fc= 1 and, thus, the amount of tissue being resorbed is only controlled in this model by the number of active BMUs, through equation (1). This is not actually true since the resorption area of a BMU may vary, likely under the influence of mechanical loading (van Oers et al., 2008b) and thus fccould differ from 1, but, the variation of the number of active BMUs has the same effect and is controlled by mechanical loading as well.fbis defined as a function of the mechanical stimulus: greater than 1 in overload and lower than 1 in disuse (see (Mart´ınez-Reina et al., 2009) for further details). It is then assumed that more osteoblasts will be differentiated with an increasing stimulus (Martin et al., 1998). The factor (1 −p) is introduced in (6a) to account for the fact that BMUs are resorbing a porous tissue, with porosity p. Finally, ARand AFare two variables having dimensions of area. They are introduced to measure the amount of tissue being resorbed and formed at a certain stage of the resorption and formation processes, respectively. In the old model ARand AFwere constants and simply equal to the cross-sectional area of the tunnel resorbed by the BMU and the ring-like cylinder formed afterwards, respectively. The integrands in (6) give the contribution of those BMUs activated between time t′and time t′+dt′, contribution which is integrated over the recent history, thus taking into account only those BMUs which are still active. This is automatically done in (6) since AR(t”) and AF(t”) are identically zero for t” over a certain value. More precisely, they become zero when the BMU finishes its activity, as will be explained later on. Any difference between the volumes of tissue resorbed and formed causes a change in bone porosity. Apart from this, the activity of BMUs influences on other bone properties in several ways. For example, if that porosity is oriented, as it uses to be in cortical bone, remodelling can also modify bone anisotropy. The mineral content is also changed. More precisely, it decreases with the activity of BMUs, as they resorb old, mineralized bone, and form osteoid, which is initially not mineralized. Microstructural damage is reduced by bone remodelling as well, since the old tissue is probably damaged and the replacing tissue is initially intact. All these effects influence on the stiffness of the tissue and are accounted in the remodelling model, whose algorithm is schematically shown in figure 1. For a detailed explanation of these effects see (Mart´ınez-Reina et al., 2009). The activation of BMUs in the old model involves some simplifications 6 which are still present in the new models (both strain-based model and strain and damage-based model) and are commented next: a) Any point of cortical bone is a potential site for activation of BMUs, provided that it is placed on the bone matrix surface, mainly, vascular canals. However, in the finite element implementation of this model it was assumed that BMUs can be activated only at discrete points, the Gauss points of the mesh (Garc´ıa-Aznar et al., 2005) and regardless of being on the bone matrix surface or not. This last requirement had to be partially relaxed, since the actual bone matrix/pores interface is not modelled. Instead, porosity is homogenized and treated in the continuum mechanics sense. Then, the BMUs are activated at the Gauss points, though the model cannot distinguish whether those Gauss points are on the matrix/pore interface or not. Nevertheless, the amount of bone matrix surface is important in the activation and is accounted for through the bone specific surface, Sv, in equation (1), which is just a function of porosity. b) The BMU activation frequency is considered here as a continuous function of time and space. In this model, fractions of BMUs are continuously being activated every day in every Gauss point of the FE mesh. This way, if a fraction of BMUs is activated at a certain Gauss point in a certain day, another fraction can be activated the next day at the same Gauss point or in an adjacent one. This seems not realistic, since further activation of BMUs must be prevented in a certain volume of bone if a BMU is already progressing nearby. In such a case, it is more likely that osteoclasts were recruited to the already active BMU instead of forming a new one. A more realistic approach would have been to consider the BMU activation frequency a binary function of time and space: 0 (no BMU is activated at a certain time and position), 1 (one BMU is activated in that point at that instant). In such case, no more BMUs should be activated in the vicinity of that Gauss point until a certain time has passed. In this model the fractional approach was chosen for simplicity. These simplifications are justified since this is not a microscopical model and thus the activation, movement and remodelling effect of BMUs make sense only if averaged in a certain reference volume element (RVE), on the assumptions of continuum mechanics. 2.2. The new models. Changes introduced in the old model The new models are anisotropic like the old model. This means that the porosity and damage are treated as tensorial variables (see (Mart´ınez- 7 Reina et al., 2009) for further details). The first novelty introduced in the new models (both strain-based model and strain and damage-based model is considering the diffusion of BMUs. In the old model BMUs were activated at the Gauss points of the FE mesh, as said before, and were restricted to stay within the volume assigned to each Gauss point. The progression’s direction was fictitious since the BMU could not actually move. This direction just controlled the orientation the osteon would have if it could move, which in turn affected the local anisotropy of the tissue. In contrast to that simplistic approach, the new models do consider the movement of the BMU through the FE mesh (a trajectory can be derived) and then, the remodelling effect of the BMU is not local to the activation site any longer. 2.2.1. Resorption and formation fronts Different phases of the activity of the BMUs need to be defined (see figure 2). First, osteoclasts are being recruited at a certain location of bone. They are resorbing bone radially (rdirection, represented in phase 1 of figure 2) and progressing longitudinally (zdirection), while other osteoclasts are being recruited. This way, a resorption front, represented in dark grey, is being formed during phase 1. This phase ends when the final diameter of the tunnel, dO, is reached at the initial section. The resorption front is progressing longitudinally until the end of phase 4, when the last recruited osteoclasts undergo apoptosis. During phase 5 the rest of osteoclasts undergo apoptosis until resorption ends. On the other hand, TIdays after the resorption front was fully formed, at the end of phase 2, osteoblasts begin to secrete osteoid at the initial section. These cells deposit layers of osteoid (in light grey) creating a formation front during phase 3. This formation front is developing radially (until the haversian canal is fully formed at the end of phase 3) and progressing longitudinally during phases 4 through 6, until they reach the final section of the osteon at the end of phase 6. During phase 7 formation is receding, while the cavity created by osteoclasts is refilled with osteoid, a process which is finished at the end of phase 7. Recalling equations (6), the quantities AR(t”) · VBMU and AF(t”) · VBMU give, respectively, the volume of tissue resorbed and formed per unit time by a single BMU. Assuming that the resorption and formation fronts have a conical shape, those volumes are easily derived from the cones and truncated cones represented in figure 2. AR(t”) · VBMU ·dt” provides the volume of the cone representing the resorption front, in dark grey, while AF(t”) · VBMU ·dt” gives the volume of the truncated cone representing the formation front, in 8 light grey. Then, ARand AFare calculated as follows. AR(t”) =                  π d2 O 4(t” TR)2 Phase 1 π d2 O 4Phases 2-4 π d2 O 4[1−(t”−σL TR)2]Phase 5 (7) AF(t”) =                  π d2 O 4(1−f2 1(t”))Phase 3 π(d2 O−d2 H) 4Phases 4-6 π d2 O 4(f2 2(t”) −d2 H d2 O)Phase 7 (8) where the functions f1and f2are given by: f1(t”) = 1 −(1−dH dO)t”−TI−TR TF (9a) f2(t”) = 1 −(1−dH dO)t”−TI−TR−σL TF (9b) with dO, the diameter of the osteon and dHthe diameter of the haversian canal. Table 1 shows the values of those parameters related with the activity of BMUs which are used in the model. 2.2.2. Steering of BMUs From now on, we are only referring to cortical bone, for which the model presented in this work is aimed. Although BMUs are activated at the bone matrix surface (Frost, 1964a; Parfitt et al., 1983), at either the haversian or the volkmann canals, cortical BMUs tunnel their way through the existing bone and are not confined to progress on the bone matrix surface, as BMUs of trabecular bone are, due to the limited dimensions of trabeculi. This way, cortical BMUs may steer within the cortical layer towards certain sites where their activity is called for, in a targeted guidance (Martin, 2007). As stated before, this guidance might be driven by osteocytes undergoing apoptosis, either by disuse or by damage. Thus, damage and disuse would 9 (activation) and ending points (apoptosis) of the BMU. The histogram of figure 8 shows the distribution of this angle in both cases (the strain and damage-based model and the strain-based model) for the BMUs activated at the 256 selected points. The mean value of αin the 256 points chosen ±one standard deviation is shown in table 2 for each half and each model. The Kolmogorov-Smirnov test was used to check the normality of the samples. This test shows no significant deviation from normality in the sample taken from the distal half simulated with the strain-based model (p= .598) but a significant deviation (p < .001) in every other sample (proximal half with strain-based model and both halves with the strain and damagebased model). Thus, non-parametric tests were used to compare the samples. The U-Mann-Whitney test revealed a significant difference in the orientation of osteons between the distal and the proximal halves in both cases: strain and damage-based model (Z= 3.85, p < .001 and medium effect size r=.24) and strain-based model (Z= 4.60, p < .001 and medium effect size r=.29). The Wilcoxon signed-rank test for related samples revealed that there is no significant difference between the orientation of osteons obtained with the strain and damage-based model and with the strain-based model for the distal half (Z= 0.697, p=.486 and a small effect size r=.04) while there is a difference for the proximal half, though not significant at the level .05 (Z= 1.658, p=.097 and a small effect size r=.10). 4. Discussion Figure 5 shows the effect of the strain and damage-based model, against a variation of it, the strain-based model, similar to the old model. The trajectories of BMUs are basically the same except for the area with a high damage level (see figure 6), where there is a noticeable difference. In figure 5a, using the strain and damage-based model,Kis defined by equation (20) and thus depends on the signal level, and consequently on the damage level (recall equation (2)). In this case, BMUs are attracted and guided by damage and then, those BMUs progressing near the area with a high damage level, are steered to the periosteum where the damage is higher, repairing the damage in their way. In figure 5b, using the strain-based model,Kis forced to be null and BMUs are not attracted by damage and just follow the direction of the prevalent local stress directions as hypothesized by Burger et al. (2003). The effect of this is that BMUs repair damage more efficiently using the 16 strain and damage-based model as can be seen in figure 6 and as suggested by Martin (2007). Another undesired effect of the hypothesis of BMUs following just the prevalent strain directions (done in the strain-based model and the old model) can be seen in figure 6 (right), more precisely in the BMUs progressing near the area with a high damage level. Damage affects the stiffness tensor of the material (orthotropic in this model) and this affects the principal strain directions. The result of a damage distribution like that is an irregular orientation of the prevalent strain direction, which makes those BMUs to be somewhat disoriented, progressing in directions which have no clear purpose. It must be said, though, that this drawback derives from a numerical problem and it is not directly attributable to the hypothesis of Burger et al. (2003), but to the difficulty of implementing such hypothesis in a FE analysis. In the areas of the femur where the damage level is moderate to low, no differences can be seen in the trajectories of the BMUs. The statistical analysis confirmed this conclusion. In the distal half of the FE model, no significant difference between the orientation of the osteons obtained with each model was found, while in the proximal half (containing the area of high damage) there is a difference, though not significant. This different behaviour of the proximal and distal halves of the model was also statistically confirmed. Significant differences were found for both models. In the case of the strain and damage-based model this difference arises from the deviation of osteons to the damaged area in the proximal half, while using the strain-based model that deviation is due to the undesired effect of those “disoriented” BMUs, mentioned right above. The inclination of osteons obtained in the simulations varies between 10◦ and 30◦, being the range 10 ÷20◦the most numerous. This result agrees with other results found in the literature. Heˇrt et al. (1994) measured the inclinations of osteons in human femora, obtaining the range 5 ÷15◦. That difference may arise from the fact that these authors measured that inclination of osteons in plane sections cut in the femur. This way, they could only appreciate the plane section of the haversian canal and the angle measured is that of the projection of the osteon on the plane section, which is smaller than the actual one. To illustrate this, let us consider a bone with the long axis in the Zdirection of a certain coordinate system (see figure 9). For an osteon with a direction v, the actual angle αmust be measured in the plane πcontaining Zand the osteon. If the plane section, π′, cut to measure the inclination of the osteon did not contain the osteon, but were 17 parallel to the plane XZ, for example, the angle measured for that osteon would be α′. This angle depends on the angle β, formed by planes πand π′. If we assume that βfollows a uniform probability distribution, it can be shown that the expected value of the measured angle α′for an average angle α= 20◦would be α′= 14.3◦. Likewise, for an actual average α= 10◦, the measured average would be α′= 7.1◦. Then, the range 10◦÷20◦obtained in this work as the most likely, is transformed into the range 7.1◦÷14.3◦, which, very interestingly, is quite close to the range 5◦÷15◦measured by Heˇrt et al. (1994). Figure 10 shows some of the sections provided by Heˇrt et al. (1994). After staining with India-ink, the haversian canals are shown in black. A wide scatter can be noticed in the length of the vascular canals. Many of them are long, but some others are of medium length and some are even very short. Such a great variation in the length of the haversian canals has not been reported, what would suggest that many of the canals appearing in the picture are actually plane sections of longer canals. These canals would have been cut by a plane like π′and not by a plane containing the canal itself, needed to compute the actual inclination. This would manifest two important conclusions: firstly, that the average values of inclination of osteons provided by Heˇrt et al. (1994) would be underestimated. Secondly, it would also confirm the wide scatter of osteons’ inclination, especially in the distal part of the diaphysis (compare with table 2). This fact could be behind the differences between the results obtained here (30.9±21.5◦for the distal half of the diaphysis and 20.1±6.6◦for the proximal half) and the averages reported by Heˇrt et al. (1994) (8.1◦in the lateral side and 10.3◦in the medial side). Another interesting observation of Heˇrt et al. (1994) that can be seen in figure 11 is the presence of a typical fan distribution of vascular canals near the periosteum. This is in accordance with the results obtained with the strain and damage-based model, more precisely, with the trajectories of BMUs near the damaged area (see figure 5). These BMUs turn to the periosteum, where the damage is higher due to the bending stresses. On the contrary, the strain-based model does not reproduce this result for the BMUs passing close to the periosteum. Instead, they turn in an irregular way and are finally driven away from the periosteum, as stated before. The damage accumulation rate has a major influence on the trajectories predicted by the strain and damage-based model and is highly dependent on two things: the mechanical properties of bone and the stress state the femur 18 is subjected to. Validated experimental results dealing with bone fatigue are very scarce. In fact, only a few models can be found in the literature (one of them, the model used here) and their validity must be considered with care. Regarding the distribution of stresses in the femur it must be admitted that the loads used in this model, taken from the literature as well, are only approximate. For example, not all the muscular groups, but only the most important one, were included and besides, the boundary conditions were very simplistic. Therefore, the results obtained in this study must be analysed from a qualitative point of view, rather than quantitative. Some comments about how this model counts the activated BMUs are pertinent at this moment. Activation of BMUs is a process by which a discrete number of BMUs can be activated within a somewhat continuous region of bone (BMUs are only constrained to be activated on the bone free surface, what in cortical bone means on the surface of vascular canals). It seems plausible that the site of BMU activation should be the point with the highest value of activating signal in a certain reference volume, whichever that signal is. Once the preosteoclasts have gathered around this attracting point to form multinucleated osteoclasts, the activation of further BMUs might be presumably inhibited in the surrounding volume. Opposite to this idea, this model assumes that a continuous number of BMUs can be activated in a discrete number of points throughout the bone (the Gauss points of the finite element mesh). In addition, no matter if a fraction of BMUs was activated a certain day, the next day, another fraction of BMUs can be activated in the same point. This is a strong simplification of the problem, but, simulating the real process would require the use of a discrete activation frequency and even a microscopical model, in order to account for the effect of a single BMU in the mechanical properties of the surrounding bone and this is out of the scope of this work. Nevertheless, the simplification proposed here to circumvent this difficulty can be interpreted in the following way: if the activation of BMUs is kept high for a certain period of time in a certain Gauss point, the fractions of BMUs activated daily will have, approximately, the same effect than a discrete number of BMUs, provided that this effect is averaged over time and space. So, the simplification made in this study can be justified based on the assumptions of continuum mechanics. Therefore, the BMUs in figure 5 and in the histogram are not really whole BMUs but fractions thereof. So, figure 5 is actually representing the tendency of the BMUs progressing by these areas. In any case, if a remodelling equilibrium situation were to be modelled, in which no changes of mechanical stimulus 19 and properties of the bone take place, no significant change of that tendency of the BMUs would be encountered. Thus, the counting of BMUs done here would be accurate enough. The mechanical stimulus defined in this model (see equation (4)) is proportional to the number of load cycles. However, it has been suggested the existence of a threshold in the number of cycles for the saturation of the stimulus (Adams et al., 1997). This saturation may be due to the fact that osteocyte stimulation is much reduced after the first cycles of repetitive loading (Srinivasan et al., 2002). The effect of such saturation on the bone remodelling model proposed in this paper is to yet be analysed. Finally, it must be noted that the high number of parameters and the complexity of the bone remodelling models proposed in this paper makes the results extremely difficult to validate quantitatively. Moreover, the value of some parameters lacks validation and were chosen based on previous bone remodelling models. Others, like BMU lifespan, BMU progression rate and the length of the different phases, etc., were taken as deterministic constant values. However, some works have suggested that bone remodelling might be a stochastic process (Christen et al., 2013; Dunlop et al., 2009; Hartmann et al., 2011; Weinkamer et al., 2004) instead of a deterministic one. In addition, other works suggest that the BMU lifespan, BMU progression rate and the length of the phases might be controlled by certain biochemical signals (Buenzli et al., 2011; Pivonka and Komarova, 2012), which may make those parameters variable. Neither the randomness, nor even the variation of those parameters was considered in this work. 5. Conclusions This paper presents an anisotropic remodelling model for cortical bone based on the directional activity of BMUs. This model is an extension of a previous one (Mart´ınez-Reina et al., 2009), in which several improvements have been included in order to better account for the effect of BMUs in the remodelling activity, specially in microdamage repairing. In this regard, the hypothesis that BMUs follow the direction of the prevalent local stress (Burger et al., 2003) has been revised according to the observations of Martin (2007). This author argued that BMUs must be somehow driven by damage, given that the efficiency of damage repairing by bone remodelling observed in histomorphometric analyses was higher than expected. This higher efficiency would be explained, according to Martin (2007), by the possibility that 20 BMUs could turn along their way, attracted by certain signalling, expressed by apoptotic osteocytes present near microcracks. This coupled steering of BMUs by damage and prevalent strain direction has been included in the model proposed here. This model has been tested on a FE model of a human femur subjected to normal walking loads. The effect of assuming that BMUs are guided only by the prevalent strain directions or also by damage (as in the proposed model) has been analysed. If both variables were responsible for the guidance of BMUs, damage repairing would be more efficient and BMUs would not follow an irregular course in high damaged areas as would be the case if BMUs were only guided in the direction of the prevalent local strain. In some cases, this irregular course even makes BMUs to run away from the damaged area, what would be in contradiction with one of the main purposes of bone remodelling: microstructural damage repairing. The assumption of BMUs following the direction of the prevalent local strain works reasonably well and may explain, up to a certain point, the orientation of osteons, but fails to explain some specific phenomena like the high efficiency of BMUs in damage repairing just discussed and the fan distribution of osteons near the periosteum (Heˇrt et al., 1994), where the amount of microstructural damage can be high due to the bending stresses. Therefore, it can be concluded that it is crucial to consider the effect of damage in the steering of BMUs in a remodelling model of cortical bone. 6. Appendix A The density of apoptotic cells can be determined by the concentration of two molecules: Bax and Bcl-2 (Verborgt et al., 2002). Bax is one of the principal effector molecules commonly expressed in cells undergoing apoptosis, including osteocytes. Bcl-2 is a pro-survival molecule having a molecular structure similar to Bax but an opposite function. Bcl-2 counteracts the effect of Bax acting as a shield against macrophages. In the vicinity of a crack the concentration of Bax and Bcl-2 was reported by Verborgt et al. (2002). The variation of the net concentration, i.e., the difference between the concentration of Bax and Bcl-2, is normalized by dividing by its maximum value (found at the crack tip) to obtain the variable b: b=[Bax]−[Bcl −2] ([Bax]−[Bcl −2])max (18) 21 The variation of this normalized concentration, b, with the distance from the crack is schematized in (Martin, 2007) (see figure 12). Following this author, this net concentration acts as a signal to attract macrophages, osteoclasts in this particular case. Following Burger et al. (2003), we have assumed that osteocytes placed in a region of bone with a low mechanical stimulus will have a similar behaviour to those in the vicinity of a crack. Consequently, the attracting signal,s, which accounts for damage and disuse (see equation (3)), is assumed to be somewhat related to the net [Bax] - [Bcl-2] signal, and consequently to the density of apoptotic osteocytes. The attracting signal factor, K, was defined to vary with sin the strain and damage-based model.Kcan be thought as a measure of the normalized concentration of [Bax] - [Bcl-2], b. In a totally damaged (or completely disused) area of bone, s= 1, by definition of the attracting signal, and, logically, the concentration of [Bax] - [Bcl-2] should be at its maximum, b= 1. In this situation BMUs are assumed to be steered to the damaged/disused area regardless of the direction of the prevalent strain, that is K= 1. So, the identification of Kand bis consistent at this point (s= 1, K =b= 1). Such identification has also been made at the point (s0, K0=b0), which represents a reference or average situation. These values are obtained next. The reference situation is represented by a certain density of microcracks: a single microcrack of average length (76µm following Martin (2007)) centered in a sphere of radius r0(see figure 12). If cracks are the source of the concentration of [Bax] - [Bcl-2] and there are no other cracks within that reference volume, the distribution of bwill be radial like that represented in figure 12 and the average concentration in this volume can be calculated by integration. b0=∫V0 b(r)dV ∫V0 dV =1 4(r1 r0)3 (1 + 3 b1)+ +b1 4 (r0−r1)(r0−4r3 1 r2 0 +3r4 1 r3 0) (19) where the following values of the parameters: r0= 4 mm,r1= 1.5mm, b1= 0.2 are taken from (Martin, 2007), resulting in b0= 0.08897. Now, it can be thought, on the assumptions of continuum mechanics, that a piece of bone with a crack density equal to that of the reference situation 22 (1 crack/(4πr3 0/3) = 3.73 ·10−3cracks/mm3) has an average concentration of [Bax] - [Bcl-2] equal to b0= 0.08897 and then, through identification, K0= 0.08897. On the other hand, that crack density produces a stiffness loss of approximately 20% in canine femora under bending (Pidaparti, 2000; Burr, D.B., personal communication) and that loss of stiffness corresponds to diso = 0.0352 (see Appendix B). Substituting this value into equation (2) for a situation of normal load (ξ=ξ∗), the value s0= 0.535 is obtained with the parameters aand c=ξ∗of table 1. The point (s0= 0.535, K0=b0= 0.08897) tries to represent a normal situation, over which the attractive effect rises up considerably with damage or disuse. This rise has been assumed to occur in a linear manner, up to the point (s= 1, K = 1). Below the reference point (s0, K0=b0), the attracting effect is assumed to decay exponentially to (s= 0, K = 0). In this last situation, no attracting signal exists and the concentration of [Bax] - [Bcl-2] must be zero, that is K=b= 0. The final function K=K(s) is: K(s) =    A(eB s −1) if s < s0 K0+1−K0 1−s0 (s−s0) if s0≤s < 1(20) where Aand Bmust be chosen to ensure C1continuity at s=s0. 7. Appendix B The old model considers the directionality of damage and follows the Theory of Anisotropic Continuum Damage proposed by Cordebois and Sidoroff (1982). There, the microstructural damage is defined as a tensorial magnitude. For example, a damage tensor like the following: d=  0 0 0 0 0 0 0 0 d (21) would be produced by the accumulation of damage in an initially undamaged material subjected to uniaxial load in zdirection. The compliance tensor of the damaged material, S, is related to that of the undamaged material, ˆ S, by: S= (I−D)−1(I−D)−1ˆ S(I−D)−1(I−D)−1(22) 23 in order to maintain the symmetry (see (Cordebois and Sidoroff, 1982)). The damage accumulated by the specimens tested by Pidaparti et al. (2000), subjected to uniaxial bending, would be modelled by a damage tensor like that of equation (21). In that case, the loss of stiffness in z direction, or, alternatively, the remaining stiffness could be easily derived from (22): Ez ˆ Ez = (1 −d)2(23) where Ezand ˆ Ezare, respectively, the Young’s moduli of the damaged and undamaged material in the axial direction of the bending test. The parameter diso in equation (2) is related to the microstructural damage present in the piece of bone and is defined in (Mart´ınez-Reina et al., 2009) as diso =tr(d)/3, with tr(·) designating the trace of a tensor. Then, in uniaxial damage diso = d/3. Finally, for a loss of stiffness of 20% (or a remaining stiffness of Ez/ˆ Ez= 0.8) diso = 0.0352. 8. Acknowledgements The authors gratefully acknowledge the financial support given by the Consejer´ıa de Econom´ıa, Innovaci´on, Ciencia y Empleo of the Junta de Andaluc´ıa for the development of the project P09-TEP-5195. 9. References [1] Adams, D.J., Spirt, A.A., Brown, T.D., Fritton, S.P., Rubin, C.T., Brand, R.A., 1997. Testing the daily stress stimulus theory of bone adaptation with natural and experimentally controlled strain histories. J. Biomech. 30, 671-678. [2] Aluru, N.R., Li, G., 2001. Finite cloud method: A true meshless technique based on fixed reproducing kernel approximation. Int. J. Numer. Meth. Engng. 50, 2373-2410. [3] Beaupr´e, G.S., Orr, T.E., Carter, D.R., 1990. An approach for timedependent bone modeling and remodeling. Application: a preliminary remodeling simulation. J. Orthop. Res., 8, 662-670. [4] Bergmann, G., Graichen, F., Rohlmann, A., 1993. Hip joint loading during walking and running, measured in two patients. J. Biomech. 26, 969- 999. 24 [5] Buenzli, P.R., Pivonka, P., Smith, D.W., 2011. Spatio-temporal structure of cell distribution in cortical Bone Multicellular Units: A mathematical model. Bone 48, 918-926. [6] Burger, E.H., Klein-Nulend, J., Smit, T.H., 2003. Strain-derived canalicular fluid flow regulates osteoclast activity in a remodelling osteon-a proposal. J. Biomech. 36, 1453-1459. [7] Christen, P., Ito, K., dos Santos, A.A., M¨uller, R., van Rietbergen, B., 2013. Validation of a bone loading estimation algorithm for patient-specific bone remodelling simulations. J. Biomech. 46, 941-948. [8] Cordebois, J.P., Sidoroff, F., 1982. Damage induced elastic anisotropy, in: Boehler, J.P. (Ed.) Mechanical behavior of anisotropic solids , Proc. EUROMECH Colloquium 115. Martinus Nijhoff Publishers, The Hague, Netherlands, pp. 761-774. [9] Cowin, S.C., Hegedus, D.H., 1976. Bone remodeling, I: theory of adaptive elasticity. J. Biomech. 6, 313-326. [10] Doblar´e, M., Garc´ıa, J.M., 2001. Application of an anisotropic boneremodelling model based on a damage-repair theory to the analysis of the proximal femur before and after total hip replacement. J. Biomech. 34, 1157-1170. [11] Duarte, C.A., Oden, J.T., 1996. An h-p adaptive method using clouds. Comput. Method. Appl. M. 139, 237-262. [12] Dunlop, J.W.C., Hartmann, M.A., Br´echet, Y.J., Fratzl, P., Weinkamer, R., 2009. New suggestions for the mechanical control of bone remodeling. Calcif. Tissue Int. 85, 45-54. [13] Frost, H.M., 1964a. Dynamics of bone remodelling, in: Frost, H.M. (Ed.) Bone biodynamics. Little, Brown and Co., Boston, pp. 315-333 [14] Frost, H.M., 1964b. Laws of Bone Structure. Charles C. Thomas, Springfield. [15] Garc´ıa-Aznar, J.M., Rueberg, T., Doblar´e, M., 2005. A bone remodelling model coupling microdamage growth and repair by 3D BMU activity. Biomech. Model. Mechanobiol. 4, 147-167. 25 Phase 1: Development of the resorption front t”≤TR vBMU σL vBMU dt”ZR dO r z Phase 2: Only resorption TR< t”≤TR+TI αR Phase 3: Development of the formation front TR+TI< t”≤ ≤TR+TI+TFvBMU (TI+TR) ZF dH Phase 4: Resorption + formation TR+TI+TF< t”≤σL αF dH Phase 5: Resorption recedes σL< t”≤TR+σLZF ZR Phase 6: Resorption ended. Inversion at the final section TR+σL< t”≤ ≤TI+TR+σL Phase 7: Formation recedes TR+TI+σL< t”≤ ≤TI+TR+TF+σL Figure 2: 32 K s 1.0 K0 s01.0 Figure 3: Cut boundaries R F Figure 4: 33 A P L MA P L M a b Figure 5: a A P MA P M b Figure 6: 34 ˜ d 0 30 60 90 120 150 180 210 240 270 300 0.022 0.026 0.030 0.034 0.038 K= 0 Kof eq. (13) t(days) Figure 7: 0-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80 80-90 0 20 40 60 80 100 120 140 160 180 # BMUs (degrees) strain and damage-based model strain-based model Figure 8: 35 a bp p’ a’ v Z Y X Figure 9: Figure 10: 36 Figure 11: b r b1 1.0 r1r0 Figure 12: Figure A.1. 37 Parameter Nominal value c=ξ∗Stimulus activation parameter 0.0015 aDamage activation exponent 20 fbio Maximum activation frequency 0.05 #BMUs mm2day TRResorption period 24 days TIInversion period 8 days TFFormation period 64 days VBM U BMU progression rate 40 µm/day σLBMU’s lifespan 100 days dOOsteon diameter 0.152 mm dHHaversian canal’s diameter 0.029 mm Table 1: Distal Proximal strain and damage-based model 30.9±21.5◦20.1±6.6◦ strain-based model 29.9±19.4◦18.8±3.4◦ Table 2: 38