**1. Introduction**

Bone is a living tissue that can adapt its apparent density and internal microstructure (through the process called bone remodelling), and its shape and external dimensions (through bone modelling) as a response to different mechanical and biological stimuli. Regarding the former, it has been hypothesized that one of the goals of bone remodelling is to maintain bone as an optimal structure that supports the loads with the minimum weight [1]. Thus, bone density, and consequently, stiffness, are high in overloaded regions and low in regions with a low stress level. In other words, there is a direct relationship between density and stresses. Many bone remodelling models (BRM) have been proposed in the literature to quantify this relationship [1–5].

These BRM have been very often used to estimate the density distribution in bones from the loads they are subjected to, mainly in the human femur [6,7], but also in other bones such as the mandible [8]. This problem, that we will name here *Density Prediction Through Bone Remodelling* (DPTBR), is usually approached through the following iterative process:


**Citation:** Calvo-Gallego, J.L.; Gutiérrez-Millán, F.; Ojeda, J.; Pérez, M.Á.; Martínez-Reina, J. The Correlation between Bone Density and Mechanical Variables in Bone Remodelling Models: Insights from a Case Study Corresponding to the Femur of a Healthy Adult. *Mathematics* **2022**, *10*, 3367. https:// doi.org/10.3390/math10183367

Academic Editor: Mauro Malvè

Received: 21 August 2022 Accepted: 14 September 2022 Published: 16 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).


This iterative process is repeated until convergence of the density distribution is achieved, which usually resembles the real one with good accuracy. The BRM can be applied following two approaches, as mentioned above in step 3. In an evolutionary model,the stress/strain state determines the local change in apparent density, *ρ*˙:

$$
\rho = f(S) \qquad \text{Evolutionary model} \tag{1}
$$

where *S* is a certain mechanical stimulus related to the stress/strain state. This approach is based on the Mechanostat Theory [10], which is behind most of the BRM. This theory hypothesizes that bone adapts itself to overloads by increasing its apparent density and adapts to disuse states by decreasing it. Overload and disuse are defined in the Mechanostat Theory by certain strain ranges. The theory also establishes the existence of a so-called "lazy zone", a strain range between disuse and overload, for which no evident change in apparent density is observed or, at least, it is not significant. Implementing this approach in DPTBR problems has a major drawback, as the uniqueness of the solution is not guaranteed. This path dependence occurs due to the implementation of the lazy zone [11], and the final density distribution depends on the initially assumed distribution.

Recently, we have proposed a non-evolutionary strategy to solve DPTBR problems [9]. Instead of calculating the rate of change in density as a function of the stimulus, we calculated a priori a relationship between the stimulus *S* and the density achieved at the equilibrium (see Equation (2)) and used that relationship to assign the density to an element as a function of the local stimulus. Since this stimulus can, in turn, change with density (through the stiffness), an iterative process is still required, but the number of iterations needed to achieve convergence is much lower. Notwithstanding, the uniqueness of the solution is not guaranteed in this case. The only advantage of this approach is its higher speed of convergence.

$$
\rho = \g(S) \qquad \text{Non-evolutionary model.} \tag{2}
$$

On the other hand, the non-evolutionary approach is not suitable to predict the changes in density in a bone subjected to changes in its biomechanical environment. The evolutionary approach is preferable in this case, as it allows a real-time simulation of those changes. However, to this end, it is very important to start from a realistic initial density distribution, in order to avoid the path dependence of the solution previously mentioned.

The stimulus *S* is the variable that drives bone remodelling and must comprise biological and mechanical factors. Focusing on the latter, the amount of damage has been used as the stimulus (or a part of it) in targeted bone remodelling models [12–15]. This is based on the hypothesis that one goal of bone remodelling is to repair the microstructural damage accumulated in the bone matrix by daily activity. Other models simply apply the Mechanostat Theory to account for disuse and overload and the influence of these states on bone adaptation. To this end, these models have used a stimulus that measures the intensity of the loads, with the strain energy density (SED) being the most commonly used variable to account for that intensity (see [1–5], among many others). However, in the original Mechanostat Theory, the disuse, lazy zone and overload states were defined in terms of strain; thus, establishing the hypothesis that strain is the magnitude driving the bone response. In such case, after a change in load that could alter the level of strain, the bone microstructure (and consequently, stiffness) must be regulated to return to the homeostatic situation [12].

In this work, we will focus on the mechanical part of the stimulus—the main objective being to study which is the best variable among a series of candidates to account for the mechanical feedback in bone remodelling models. We will discuss which is the best mechanical stimulus, or equivalently, which is the best predictor of bone density, either for an evolutionary or a non-evolutionary approach. For this purpose, we have analysed the stress and strain distributions in a human femur using a FE model and the loads extracted from a gait analysis performed on the same individual, from whom we have also estimated the bone density distribution from a CT scan. Several mechanical variables calculated from the stress and strain tensors throughout the gait cycle have been proposed as candidates for the mechanical stimulus (or predictors) and will be correlated with the density distribution.

This paper is organised as follows: In Section 2.1, we provide a description of the gait analysis performed to estimate the loads acting on the femur during walking. In Section 2.2, we describe the image analysis of CT scans of a human femur, which were performed to obtain its bone density distribution. In Section 2.3, we briefly describe the methodology to estimate the stiffness tensor at every point of the femur based on the density distribution. In that section, we focus on the case that considers bone as an isotropic material, while through Sections 2.4–2.6, we develop the methodology to estimate that stiffness tensor when bone is considered anisotropic. In Section 2.7, we briefly describe how the gait cycle was simulated with a FE model. The correlation coefficients used to evaluate the relationship between bone density and the variables used as candidates for mechanical stimuli (or predictors) are defined in Section 2.8. In Section 2.9, we define those predictors. In Section 3, we provide the correlations between the bone apparent density and the proposed predictors. We discuss the results of these correlations in Section 4 and highlight the conclusions of this study in Section 5. Finally, the Appendices contain a more detailed description of the procedures used in the Section 2.

### **2. Materials and Methods**

The procedure followed here is summarised in Figure 1 and is explained as follows: First, a gait analysis was performed to estimate the forces exerted by those muscles inserted in the femur and joint reactions at the knee and the hip in the subject under study. A CT scan of the subject's femur was taken and the greyscale value was related to the bone apparent density using a linear relationship. With this, a FE model of the femur was built and the loads, estimated through a gait analysis performed on the same subject, were applied to the FE model in order to obtain the distribution of stresses and strains throughout the gait cycle. A set of variables, defined in Sections 2.7 and 2.9, was assessed from the temporal evolutions of stresses and strains. Finally, the bone apparent density estimated from the CT scan was correlated with these variables.

### *2.1. Gait Analysis and Subject Data*

The gait analysis was performed at the Motion Analysis Laboratory of the Department of Mechanical Engineering and Manufacturing of the Universidad de Sevilla. A Vicon ®system of 12 infra-red cameras was used to record motion at a frequency of 100 Hz along with 2 AMTI force platforms to record the ground reaction forces at a frequency of 1000 Hz. The marker placement protocol employed was the modified Cleveland protocol [16].

The subject under study was an adult male of 27 years old, with no reported pathologies, 1.85 m tall and weighing 75 kg, who walked at a freely chosen forward speed. The participant signed an informed consent form prior to the recording of the measurements. The study protocol was approved by a medical ethics committee through the Andalusian Biomedical Research Ethics Platform (approval number 20151012181252).

**Figure 1.** The summary of the procedure followed in this work. One male subject was selected for the study (**A**). A CT scan (**D**) was used to build a FE model of the subject's femur (**F**). The gait analysis (**B**) perfomed on the same subject was used (**C**) to estimate the muscle forces (MF) and joint reaction forces (HJF: hip joint force, KJF: knee joint force) applied to the FE model, together with isostatic boundary conditions. The mechanical properties of the bone were estimated from the density, which was estimated, in turn, from the CT scan (**E**). The FE simulation of the gait cycle yielded the temporal evolution of stresses and strains (**G**) throughout the cycle. A set of variables (**I**) were derived from those temporal evolutions. Finally, these variables were correlated (**H**) to the bone apparent density in order to evaluate a good predictor of density for bone remodelling models.

The recorded experimental data were processed using OpenSim software [17]. The biomechanical model implemented was the *Gait2392*, available in the OpenSim library. This model was developed to perform 3D gait analysis and consists of 23 degrees of freedom and 92 actuators that simulate the action of muscle forces. However, for the sake of simplicity, the subtalar and metatarsophalangeal joints were blocked in this work, so that the foot was considered a single biomechanical segment. The model was scaled to fit the subject's morphology from the recorded experimental data. The mass ratio of each segmen<sup>t</sup> was assumed constant during scaling. In this model, muscle attachment points were placed where OpenSim locate them by default, using numerical approximation [18] of cadavers' data [19,20].

The forces applied to the FE model were obtained in two steps. First, the inverse dynamic problem was solved using the kinetics experimental data recorded in the laboratory. From these results, the time evolutions of the muscle forces were calculated solving a forcesharing problem through a static optimisation algorithm that considered the dynamics of muscle contraction and activation [17,21]. The results were collected for those muscles inserted in the femur and are provided in the Supplementary Materials. Additionally, joint reaction forces at the hip and knee were estimated using the algorithm proposed by Steele et al. [22]. A detailed description of this procedure is included in Appendix B.

#### *2.2. CT Scan of the Femur, FE Model and Density Distribution*

A CT scanner (LightSpeed16, GE Medical Systems, Milwaukee, WI, USA; 120 kVp, reconstructed with bone Plus kernel, 1.25 mm slice thickness) was used to estimate bone density on the right femur of the same individual on which the gait analysis was conducted. The software Abaqus (version 2020, SIMULIA, Dassault Systémes, Madrid, España) was used to run FE analysis. The 3D geometry of the femur was reconstructed from the CT scan and meshed with 4-node linear tetrahedral elements (C3D4 in Abaqus element library). The final mesh consisted in 339,168 elements and 64,757 nodes. In order to check the convergence of the FE solution, we compared the results with a different mesh, built with 10-node quadratic tetrahedral elements (C3D10). This new mesh was obtained from the previous one by simply placing mid-side nodes in the edges of the former elements, resulting—obviously—in the same number of elements and 480,480 nodes.

To assign a value of density to each element of the mesh, a linear relationship between the greyscale value and the bone apparent density was used, as in [23–26]:

$$
\rho\_{\varepsilon} = A + B \cdot \lg\_{\varepsilon} \tag{3}
$$

where *ge* and *ρe* are the greyscale value and the density estimated for element *e*, respectively. Two pairs of greyscale–density values are needed to define the linear relationship (constants *A* and *B*). These points are usually obtained through calibration of the CT scan. However, this calibration was not available in our case, and therefore, we needed to make two assumptions. Thus, as the first pair, an apparent density value of 2.1 g/cm<sup>3</sup> [4] was assigned to the maximum greyscale value of the CT scan (255), i.e., to the densest cortical bone. The second point chosen was that corresponding to a null bone apparent density, which is sometimes called a grey level (or greyscale) threshold (GLT). It can be seen in Equation (4) that *ρe* = 0 for *ge* = *GLT*. This lower limit is usually made to correspond to bone marrow, which was assumed here to be placed inside the diaphyseal canal. However, the greyscale value inside the canal varied in the range 70–90, thus, making it impossible to assign a unique and reliable value to GLT. For this reason, three cases were studied, assuming different values for GLT, namely: 70, 80 and 90. This uncertainty in the choice of GLT affects the slope of the linear relationship, which is then:

$$\rho\_{\varepsilon} = 2.1 \cdot \frac{g\_{\varepsilon} - GLT}{255 - GLT}.\tag{4}$$

The greyscale distribution was mapped from the CT scan to the FE mesh using the software *Bonemat*, and Equation (3) was used to convert the greyscale value into bone apparent density. Those elements with a greyscale value below GLT were assigned a very small density (0.001 g/cm3), thus, yielding a negligible—although, not null—stiffness, something necessary to avoid convergence problems. Figure 2 shows the distributions of density obtained for the three GLT. It can be seen that those three GLT led to distributions of the estimated density which are similar, in general, with the exception of the proximal region and the thickness of the cortical layer in the diaphysis. As GLT rises, the volume identified as marrow increases, and this makes GLT = 90 produce a thinner cortical layer and a slightly underestimated density in the proximal region. Nonetheless, it can be noted in Figure 2 that the uncertainty introduced in the density distribution by GLT is not too important, in the range of 70–90. Moreover, we will show later that the conclusions of this study are the same regardless of the choice of GLT.

The external boundary of the bone was not perfectly defined in some slices of the CT scan where the cortex was too thin, since the average greyscale value between the background and the periosteum produced a cortex of intermediate density. For this reason, the model was covered with a layer of shell elements to simulate the cortex, with a thickness of 1 mm, a Young's modulus equal to 19 GPa [27] and a Poisson's ratio *ν* = 0.32 [28].

**Figure 2.** Distribution of bone density estimated from Equation (3) for three GLT: in the whole femur (**left**), in a diaphyseal cross-section (**right**).

### *2.3. Estimation of the Stiffness Tensor*

Isotropic and anisotropic models were considered. In the former, the stiffness tensor of bone at each material point can be estimated through the density obtained from the CT scan, by using one of the numerous correlations between the elastic constants and the apparent density that can be found in the literature. For example, Jacobs proposed the following [28]:

$$E\_{\text{Jacobi}}\left(\text{MPa}\right) = \begin{cases} 2014\,\rho^{2.5} & \text{if } \rho \le 1.2 \text{ g/cm}^3\\ 1763\,\rho^{3.2} & \text{if } \rho > 1.2 \text{ g/cm}^3 \end{cases} \tag{5a}$$

$$\nu = \begin{cases} \ 0.2 & \text{if } \rho \le 1.2 \text{ g/cm}^3 \\ \ 0.32 & \text{if } \rho > 1.2 \text{ g/cm}^3 \end{cases} \text{.} \tag{5b}$$

Hernandez et al. [29] proposed another correlation for *E*, based on the ash fraction, *α* (a variable used to measure the mineral content of bone tissue), and on the bone volume fraction (or bone volume per total volume), *BV*/*TV*. If we express *BV*/*TV* = *ρρ*ˆ , with *ρ* and *ρ*ˆ being the apparent density and the tissue density, respectively, then the correlation proposed by Hernandez et al. reads:

$$E(\text{MPa}) = 84 \,\text{\AA} \,\text{\(\frac{\rho}{\hat{\rho}}\)}^{2.58} a^{2.74}.\tag{6}$$

If typical values of ash fraction *α* = 0.68 and tissue density *ρ*ˆ = 2.31 g/cm<sup>3</sup> [29] are used, Equation (6) becomes:

$$E\_{\text{Hermudez}}(\text{MPa}) = 3388 \,\rho^{2.58},\tag{7}$$

which produces an estimation of the Young's modulus up to 70% higher than Equation (5a) depending on the density (see Figure 3). Each correlation was used in the corresponding isotropic model, named, respectively, *IsoJ* (Equation (5a)) and *IsoH* (Equation (7)) after Jacobs and Hernandez. A constant Poisson's ratio *ν* = 0.3 was used in conjunction with Equation (7) in model *IsoH*.

**Figure 3.** Comparison of the correlation provided by Jacobs [28] and the one derived from the correlation of Hernandez et al. [29] for the Young's modulus as a function of bone apparent density.

Bone is actually an anisotropic material with a dependence of the mechanical properties on the direction and the type of tissue. Particularly, cortical bone is usually modelled as a transversely isotropic material [30] and trabecular bone as an orthotropic material [31]. Consequently, its anisotropy must be taken into account in the estimation of the stiffness tensor, and so was carried out in the anisotropic model, named here *AnisoH* (note that the H refers to the fact that the anisotropic model relies on Hernandez correlation, as explained later in Section 2.4).

Some bone remodelling models have been proposed to predict not only the adaptation of bone apparent density but also of its anisotropy [4,5]. The model developed by Doblaré and García [4] was also used to predict the distribution of anisotropy in the same manner as DPTBR simulations are used to predict the distribution of bone density. In fact, these authors predicted both distributions simultaneously. The starting point was an isotropic material with a uniform distribution of density. Their BRM adapted bone density and anisotropy, the latter through the fabric tensor, and both variables were used to update the stiffness tensor. Convergence was deemed when both density and anisotropy remained constant between simulations. A brief summary of this model is provided next, in Section 2.4.

Given that DPTBR simulations were shown to be path-dependent [11], we estimated the density distribution from the CT scan and used the bone remodelling model developed by Doblaré and García [4] to estimate the anisotropy. This required a slight variation of the original model, explained in Section 2.5. Besides, another variation of the original model was used to account for time-varying loads, such as walking. This variation was introduced by Ojeda [32] and is presented in Section 2.6.

The objective of comparing the three models referred to above (*IsoJ*, *IsoH* and *AnisoH*) was to rule out that the assumption made for the constitutive model is forcing a certain correlation between the density and the predictors.

### *2.4. Anisotropic Bone Remodelling Model Based on Continuum Damage Mechanics—Model by Doblaré and García*

A brief description of the anisotropic bone remodelling model developed by Doblaré and García [4] is given next. Consulting the original paper is advised for a more detailed description of the model. This model is an extension to the anisotropic case of the model developed by Beaupré et al. [3]. It applies the Mechanostat Theory by defining a stepwise

linear relationship between the mechanical stimulus, Ψ*t*, and the bone formation/resorption rate, *r*˙ (see Figure 4). The reference value of the stimulus, <sup>Ψ</sup><sup>∗</sup>*t* defines a region of width 2w called lazy zone (analogous to the "adapted-window" of the Mechanostat Theory), where no net change of bone density is produced. The value of *r*˙ determines the temporal evolution of apparent density, *ρ*˙. In turn, apparent density determines the variation of the Young's modulus through correlations such as (5a) or (7).

**Figure 4.** The Mechanostat Theory as interpreted by Beaupré et al. [3]. This law establishes a stepwise linear relationship between the bone formation/resorption rate, *r*˙, and the difference between the mechanical stimulus, Ψ*t*, and a reference value of that stimulus, <sup>Ψ</sup><sup>∗</sup>*t*.

The model developed by Doblaré and García was based on the theory of Continuum Damage Mechanics (CDM), where the stiffness tensor of the damaged material, **C**, is obtained from the tensor of the the non-damaged material, **C0**, and damage:

$$\mathbf{C} = (1 - D)\mathbf{C}\mathbf{o}\_{\prime} \tag{8}$$

where *D* is the damage variable, null for an intact material and equal to 1 for a completely damaged or failed material. In the extension of CDM to the anisotropic case introduced by Cordebois and Sideroff [33], the scalar damage *D* is replaced with a damage tensor **D** and the resulting material is orthotropic, with the principal directions of orthotropy aligning with the principal axes of the damage tensor **D**. Damage is understood as a measure of porosity and the directionality of that porosity and both are incorporated into the model jointly, by following the idea suggested by Cowin [34] for the fabric tensor, **H**. Therefore, the undamaged material is an ideal situation of a perfectly isotropic bone with null porosity. The damage and fabric tensors are related by:

$$\mathbf{D} = \mathbf{1} - \mathbf{H}^2,\tag{9}$$

with **1** being the second order identity tensor. Equation (9) leads to the following relationship between the components of the compliance tensor in the principal directions and the eigenvalues of the fabric tensor, *hi*:

$$\begin{aligned} \frac{1}{E\_i} &= \frac{1}{\hat{E}} \frac{1}{h\_i^4} \\ -\frac{\nu\_{ij}}{E\_i} &= -\frac{\hat{\upsilon}}{\hat{E}} \frac{1}{h\_i^2} \frac{1}{h\_j^2} \\ \frac{1}{G\_{ij}} &= \frac{1+\hat{\upsilon}}{\hat{E}} \frac{1}{h\_i^2} \frac{1}{h\_j^2} \end{aligned} \tag{10}$$

where *E* ˆ and *ν*ˆ are, respectively, the Young's modulus and Poisson's ratio of the bone with no porosity. These values can be obtained through correlations (5) or (7) for an apparent density *ρ* equal to the density of the bone matrix, *ρ*ˆ, which is assumed constant. In this particular model*AnisoH*, we assumed *ρ*ˆ = 2.1 g/cm<sup>3</sup> and chose the Hernandez correlation (7), together with *ν* = 0.3. The influence of porosity and directionality are factorised in this model by redefining the fabric tensor as follows:

$$\mathbf{H} = \left(\frac{\rho}{\delta}\right)^{\beta/4} A^{1/4} \,\mathbf{\hat{H}}^{1/2} \,\text{\AA} \tag{11}$$

where **H** ˆ is the normalised fabric tensor. This tensor is normalised by imposing det(**H**<sup>ˆ</sup> ) = 1 in order to account only for the directionality of the pores. Additionally, *β* is the exponent of the apparent density in the correlations (5a) or (7) and *A* is a parameter introduced to ensure that the formulation reproduces the isotropic bone remodelling model developed by Beaupré et al. [3] if it is applied to an isotropic case. *A* can here be considered a constant. As stated before, the quotient *ρ*/*ρ*<sup>ˆ</sup> in Equation (11) is equal to the bone volume fraction, *BV*/*TV* = 1 − *p*, with *p* as the porosity. The mechanical stimulus is defined in this model through the tensor, **Y**:

$$\mathbf{Y} = 2\left\{ 2\mathbf{\hat{G}}\operatorname{sym}\left[ (\mathbf{H}\boldsymbol{\varepsilon}\mathbf{H})(\mathbf{H}\boldsymbol{\varepsilon}) \right] + \boldsymbol{\lambda}\operatorname{tr}(\mathbf{H}^2\boldsymbol{\varepsilon})\operatorname{sym}(\mathbf{H}\boldsymbol{\varepsilon}) \right\},\tag{12}$$

where *G* ˆ and *λ* ˆ are the Lamé constants corresponding to the cortical bone with no porosity and tr(•) and sym(•) represent the trace and symmetric part of a tensor, respectively. Doblaré and García defined another tensor, **J**, to quantify the relative influence of the spherical and deviatoric parts of the stimulus:

$$\mathbf{J} = \frac{1-\omega}{3}\operatorname{tr}(\mathbf{Y})\mathbf{1} + \omega \operatorname{dev}(\mathbf{Y}),\tag{13}$$

where **1** is the identity tensor, dev(•) represents the deviatoric part of a tensor and the anisotropy factor, *ω*, weights the importance of the anisotropy of the stimulus in the model. This factor ranges from *ω* = 0, which means that the model only depends on the isotropic component of the stimulus, to *ω* = 1, which produces the maximum level of anisotropy. The same value used by Doblaré and García [4] was used here (*ω* = 0.1). Two functions *gr* and *gf* are proposed to establish the remodelling criteria. These functions depend on the stimulus **J** and are allowed to distinguish the formation, resorption and lazy zones, as carried out in Figure 4. For that reason, those functions also depend on the reference value of the stimulus, <sup>Ψ</sup><sup>∗</sup>*t* , and the width of the lazy zone, through w. Their expressions are quite complex and can be consulted in [4]. The remodelling criteria are given by the following conditions:

$$\begin{array}{llll} & \mathcal{g}\_{f}(\mathbf{J}, \mathbf{\varvarproj}\_{t}^{\*}, \mathbf{w}) \le 0 & \mathcal{g}\_{r}(\mathbf{J}, \mathbf{\varproj}\_{t}^{\*}, \mathbf{w}) > 0 & \text{resorption}; \\ & \mathcal{g}\_{r}(\mathbf{J}, \mathbf{\varproj}\_{t}^{\*}, \mathbf{w}) \le 0 & \mathcal{g}\_{f}(\mathbf{J}, \mathbf{\varproj}\_{t}^{\*}, \mathbf{w}) > 0 & \text{formation}; \\ & \mathcal{g}\_{f}(\mathbf{J}, \mathbf{\varproj}\_{t}^{\*}, \mathbf{w}) \le 0 & \mathcal{g}\_{r}(\mathbf{J}, \mathbf{\varproj}\_{t}^{\*}, \mathbf{w}) \le 0 & \text{lazy zone.} \end{array} \tag{14}$$

Based on the fulfilment of the corresponding criterion, the variation of the fabric tensor **H**, that accounts for the variation of anisotropy (through tensor **H**ˆ ) and porosity (see Equation (11)), is provided by:

$$\begin{aligned} \mathbf{\hat{H}} &= k\_1 \frac{\not\not\rho}{\rho} \mathbf{J}^{-3} \; \mathbf{\hat{o}} & \quad \text{resorption};\\ \mathbf{\hat{H}} &= k\_2 \frac{\not\not\rho}{\rho} \mathbf{J} \mathbf{\hat{o}} & \quad \text{formation};\\ \mathbf{\hat{H}} &= 0 & \quad \text{lazy zone}; \end{aligned} \tag{15}$$

where the tensor *ω***<sup>ˆ</sup>** is introduced to simplify the expression, as follows:

$$
\hat{\omega} = \frac{1 - 2\omega}{3} \mathbf{1} \otimes \mathbf{1} + \omega \mathbf{1} \tag{16}
$$

with **I** being the fourth-order identity tensor. The factors *k*1 and *k*2 in Equation (15) depended on several parameters in the original model. One of those parameters is *r*˙, so that the amount of formed or resorbed tissue modifies the fabric tensor through porosity.

#### *2.5. Modification of the BRM to Maintain Density Constant*

In the original model, Doblaré and García used **H ˙** to assess the variation of porosity and anisotropy, but in this work, since the density is known from the CT scan, we have forced it to remain constant and that is the reason why *k*1 and *k*2 can be assumed as constants. In such case, by deriving Equation (11):

$$\mathbf{\hat{H}} = \frac{1}{2} \left( \frac{\rho}{\rho} \right)^{\beta/4} A^{1/4} \mathbf{\hat{H}}^{-1/2} \mathbf{\hat{H}} \tag{17}$$

and given that **H** ˆ must remain normalised (det(**H**<sup>ˆ</sup> ) = 1), we finally adopted:

$$
\hat{\mathbf{H}} = c \,\hat{\mathbf{H}}^{1/2} \,\hat{\mathbf{H}},\tag{18}
$$

where *c* is a constant. This expression can be used in an Euler forward integration algorithm to yield:

$$
\hat{\mathbf{H}}(t\_{\bar{j}+1}) = \hat{\mathbf{H}}(t\_{\bar{j}}) + c \,\hat{\mathbf{H}}^{1/2}(t\_{\bar{j}}) \,\hat{\mathbf{H}}(t\_{\bar{j}}) \,\Delta t,\tag{19}
$$

where *tj*+<sup>1</sup> and *tj* are two consecutive integration steps, the time step Δ*t* = 1 is chosen and *c* is the constant necessary to enforce the condition det **Hˆ** (**tj**+**1**) = 1. The simulation started from an initially isotropic material (**Hˆ** (*<sup>t</sup>*1) equal to the identity tensor) and was stopped when the norm of the fabric tensor averaged for all the elements was almost invariable between iterations, i.e.:

$$\frac{\sum\_{\varepsilon=1}^{n} \hat{\mathbf{H}}\_{\varepsilon}(t\_{j+1}) : \hat{\mathbf{H}}\_{\varepsilon}(t\_{j+1}) - \sum\_{\varepsilon=1}^{n} \hat{\mathbf{H}}\_{\varepsilon}(t\_{j}) : \hat{\mathbf{H}}\_{\varepsilon}(t\_{j})}{\sum\_{\varepsilon=1}^{n} \hat{\mathbf{H}}\_{\varepsilon}(t\_{j}) : \hat{\mathbf{H}}\_{\varepsilon}(t\_{j})} < 0.001. \tag{20}$$

#### *2.6. Modification of the BRM to Consider Time-Varying Loads*

Doblaré and García [4] and Beaupré et al. [6] applied their models to estimate the bone density distribution in a human femur by applying the normal walking loads. These authors considered three instants of the gait cycle and treated those instants as independent loads. This procedure does not seem very plausible as they are not independent but part of the same load. For that reason, the procedure was modified by Ojeda [32] to treat the gait cycle as a single load. Moreover, the particularity of time-varying loads is taken into account with this modification.

As stated by Carter et al. [35], bone remodelling depends on the maximum stresses that the bone withstands throughout its load history. Thus, the peaks of mechanical stimulus

reached in time-varying loads would control the bone remodelling response, and it is important to note that these peaks can be reached at different instants at each bone site. Let us consider, for example, the temporal evolution of the mechanical stimulus shown in Figure 5, with three different activities. The remodelling response is assumed to depend on the maximum stimulus representative of each activity (black dots in Figure 5). The cycles are grouped by the level of stimulus and an average cycle must be chosen as representative of a certain activity. In Figure 5, A1, A2 and A3 represent, respectively, a high, medium and low intensity load. The mechanical stimulus must be obtained by superimposing the effect of all activities [35], but let us consider for a moment that only one of those loads is applied. In such case, A1 would stimulate formation, A3 resorption and A2 would produce no net change of bone mass.

**Figure 5.** Remodelling criteria with different loads. The horizontal lines limit the zones of formation (top line) and resorption (bottom line). If the peak is over the top line, then *gMf* > 0 and *gmr* < 0 and formation occurs. If the peak is below the bottom line, then *gMf* < 0 and *gmr* > 0 and resorption occurs. If the peak lies between both lines, *gMf* < 0 and *gmr* < 0 and neither formation nor resorption occurs (lazy zone).

The peaks of the stimulus coincide with the maximum of the formation criterion function, *gf* , which is proportional to the stimulus. Those peaks are termed here, *gMf* . Since the resorption criterion function, *gr*, is inversely proportional to the stimulus [4], the local minimum of this function, *gmr* , also coincides with the peaks of the stimulus (it must be noted that *gMf* and *gmr* do not necessarily coincide. They are simply reached at the same instant). Analogously, the valleys of the stimulus coincide with the minimum of *gf* and the maximum of *gr*, termed here *gmf*and *gMr* , respectively (red dots in Figure 5).

The procedure to analyse the bone remodelling process for time-varying loads begins with the calculation of the stimulus (and thus, of *gf* and *gr*) at every point of the FE mesh throughout the cycle, in order to capture the peaks *gMf* and *gmr* for each element. Based on the ideas of Carter et al. [35], Ojeda assumed that only the peaks are important from the bone remodelling perspective. Therefore, the activities plotted in solid and dashed lines in Figure 5 would lead to the same bone remodelling response, regardless of the valleys.

The criterion must identify the peaks and place them in one of the three regions: formation, resorption or lazy zone. Thus, the remodelling criteria (14) are replaced by:

$$\begin{array}{ll} \mathcal{g}\_f^M > 0 & \text{formation:}\\ \mathcal{g}\_r^m \le 0 & \mathcal{g}\_f^M \le 0 & \text{lazy zone;}\\ \mathcal{g}\_f^M \le 0 & \mathcal{g}\_r^m > 0 & \text{resorption.} \end{array} \tag{21}$$

It is important to highlight that *gMf* and *gmr* are not necessarily reached at the same instant for all the elements. This is the reason why a detailed description of the cycle is required in the modification proposed by Ojeda, as the remodelling response at each element could be driven by the stress-state reached at a different time point. For the same reason, the simplification consisting in considering three instants of the cycle as independent loads is not valid.

#### *2.7. FE Simulation of the Gait Cycle—Temporal Evolution of Stresses and Strains in the Femur*

The temporal evolution of stresses and strains in the femur during the gait cycle can be obtained by solving an elastic problem. Let **Ω** ⊂ R<sup>3</sup> be an open bounded domain and Γ = *∂***Ω** be its boundary, assumed to be Lipschitz continuous and divided into two disjoint parts Γ*D* and Γ*N* where Dirichlet and Neumann boundary conditions are applied, respectively. We denote, by **x** = *xi* **e***i*, a generic point of **Ω** and **n**(**x**) = *ni* **e***i* as the outward unit normal vector to Γ at a point **x**. The Einstein summation notation was adopted and a Cartesian basis **e***i* (*i* = 1, 2, 3) can be used without loss of generality.

Let **u** = *ui* **e***i*, *σ* = *<sup>σ</sup>ij* **e***i* ⊗ **e***j* and *ε*(**u**) = *εij* **e***i* ⊗ **e***j* denote the displacement field, the stress tensor and the linearised strain tensor, respectively. **e***i* ⊗ **e***j* (*i*, *j* = 1, 2, 3) represents the Cartesian tensorial basis. Let **b** = *bi* **e***i* denote the known vector of body forces, **t** = *ti* **e***i* the known vector of surface traction forces at Γ*N* and *δ* = *δi* **e***i* the known displacements at Γ*D*.

The Momentum Conservation Principle states:

$$
\sigma\_{\ddot{\imath},\ddot{\jmath}} + b\_{\dot{\imath}} = \rho \,\dddot{u}\_{\dot{\imath}\prime} \tag{22}
$$

where , *j* denotes the partial derivative *<sup>∂</sup>*/*∂xj* and (•) denotes the second derivative with respect to time. The right-hand side of Equation (22) represents the inertial force per unit volume. Linearised strains are related to displacements by:

¨

$$
\varepsilon\_{i\bar{j}}(\mathbf{u}) = \frac{1}{2} \left( \mu\_{i,\bar{j}} + \mu\_{j,\bar{i}} \right) \qquad i, j = 1, 2, 3 \tag{23}
$$

and finally, strains and stresses are related by the constitutive equation, which for linear elastic materials reads:

$$
\sigma\_{ij} = \mathbb{C}\_{ijkl}\,\varepsilon\_{kl} \qquad i, j, k, l = 1, 2, 3,\tag{24}
$$

where *Cijkl* are the components of the fourth-rank stiffness tensor **C**. In the case of an isotropic material, this tensor is completely defined by the Young's modulus, *E*, and the Poisson's ratio, *ν*, which are expressed as the functions of the density in the case of bone (see Equation (5)). In general anisotropic materials, this tensor has 21 independent elastic constants. In our anisotropic case, bone is assumed to be an orthotropic material, and the compliance tensor (inverse of the stiffness tensor) is given as a function of the fabric tensor (recall Equation (10)). The elastic problem is completed with the Dirichlet and Neumann boundary conditions, respectively, applied at Γ*D* and Γ*N*:

$$\mathbf{u} = \mathbf{\bar{\theta}} \quad \text{at} \quad \Gamma\_D \tag{25a}$$

$$
\sigma\_{\vec{\text{ij}}} n\_{\vec{\text{j}}} = t\_i \quad \text{at} \ \Gamma\_{\vec{\text{N}}}.\tag{25b}
$$

This boundary value problem rarely has an analytical solution, and hence, it is usually solved by means of the FEM, as carried out here.

At this point, the FE model of the femur had a complete definition of density, and consequently, of stiffness. Next, the muscle loads and joint reaction forces, estimated in the gait analysis, were applied as external forces. Muscle loads were applied as concentrated loads at the insertion points and with the direction defined by the insertion and origin points (taken from OpenSim [22]). The joint reaction forces were applied as concentrated loads on the corresponding articular surfaces, i.e., the hip reaction force on the surface of the femoral head and the knee reaction force on the surface of the epicondyles or of the epicondylar fossa. In both cases, the node of application at each instant was calculated as follows: A line was defined as passing through the corresponding joint centre (defined in OpenSim [22]) and with the direction of the reaction force at that instant. The reaction force was applied at the node closest to the intersection of the articular surface with that line. Furthermore, the minimum number of displacement boundary conditions (isostatic) was applied to restrain the rigid body motion of the FE model. The loads were varied over time, as a result from the gait analysis, but a quasi-static analysis was performed by disregarding the inertial forces in the FE simulation (*ρ u*¨*i* in Equation (22)). Nonetheless, we must note that the loads estimated with OpenSim arise from enforcing the equilibrium of all the external forces, including the inertial ones. For that reason, these inertial forces were indirectly considered in the simulations. The FE model, including the loads and the boundary conditions, is provided in the Supplementary Materials.

This pseudo-static analysis provided the temporal evolution of stresses and strains at each element *e* of the mesh (in fact, it is obtained at each integration point in full integration elements. In our case, C3D4 elements have only one integration point, which can be, therefore, identified with the element. In the case of C3D10 elements, the variables were evaluated at the centroid of the element) and at every instant *i* of the gait cycle. Several variables were derived from the stress and strain tensors (see Section 2.9). For a certain variable proposed as stimulus *S*, its value was calculated at each instant *i* and for each element, *e*, thus, yielding *Si e*. Then, the following maximum, minimum and amplitude were defined to represent its evolution throughout the gait cycle:

$$S\_{\varepsilon}^{M} = \max\_{i} (S\_{\varepsilon}^{i}) \tag{26a}$$

$$S\_{\varepsilon}^{m} = \underset{i}{\text{Min}} (S\_{\varepsilon}^{i}) \tag{26b}$$

$$\mathcal{S}\_{\varepsilon}^{A} = \mathcal{S}\_{\varepsilon}^{M} - \mathcal{S}\_{\varepsilon}^{m}. \tag{26c}$$

Note that the definition of *S<sup>M</sup> e* is related to the aforementioned hypothesis of Carter, according to which the local remodelling response would depend on the peak of stress that a bone site withstands throughout its load history [35]. As stated before, the peaks of stress do not necessarily occur at the same time for all bone sites. Therefore, considering only the loads of a single instant of the cycle is not enough to analyse the remodelling response of the whole bone. Even considering three instants of the cycle, as carried out in [4,6], is not enough. A detailed description of the gait cycle is required and the modification of the BRM presented in Section 2.6 is related to this idea. The variable, *S<sup>M</sup> e* , generalises this concept of the peak of stress to the peak of stimulus. As an alternative to the maximum of the stimulus throughout the cycle, the amplitude is considered in *S<sup>A</sup> e*.

### *2.8. Correlation Coefficients*

The Spearman and Pearson correlation coefficients were used to assess the statistical dependence between the bone density and the variables defined later in Section 2.9. Based on the definitions made in Equation (26), the following coefficients were calculated, taking each element as a point of the sample:

• R *M j* , between the maximum throughout the cycle *S<sup>M</sup> e* and the apparent density *ρ<sup>e</sup>*. • R*A j*, between the amplitude throughout the cycle *S<sup>A</sup> e* and the apparent density *ρ<sup>e</sup>*.

where *j* stands for *P* (Pearson) or *S* (Spearman). Weighted coefficients were calculated to account for the relative importance of a sample point based on the volume of the element. The expressions of the weighted correlation coefficients are given in Appendix A. The Spearman correlation coefficient is a non-parametric measure of rank correlation and it assesses how well the relationship between two variables can be described by a monotonic function. In particular, it evaluates if there is a concordance between the highest density and the highest values of a predictor variable. The Pearson coefficient is a parametric measure of correlation between two variables that assesses if they are related by a specific function. In particular, the linear, quadratic and power correlation coefficients were calculated for those variables that yielded a high Spearman coefficient, in order to confirm a strong correlation and to evaluate the best function relating the variable to bone density. Despite that strain energy density (SED) did not yield a high Spearman coefficient, it was also correlated with bone density through the Pearson coefficient and using those three functions, for reasons that will be explained later.

It is important to note that concentrated loads or displacement boundary conditions can produce spurious stress concentrations in the FE model. For this reason, the elements closest to the nodes where loads or displacement boundary conditions were applied up to a distance of two elements in all directions were removed from the correlations.

#### *2.9. Evolution of Stresses and Strains—Definition of Predictor Variables*

The set of predictor variables analysed in this work includes some variables that measure the magnitude of the stress or the strain tensor and the SED that accounts for the magnitude of both tensors in a single variable. The principal stresses are named here, *σ*1 ≥ *σ*2 ≥ *σ*3, and analogously, *ε*1 ≥ *ε*2 ≥ *ε*3. The maximum and minimum principal stresses ( *σ*1, *σ*3) and strains (*<sup>ε</sup>*1, *ε*3) are proposed as predictor variables. The maximum tensile stress is defined as: 

$$
\sigma\_t = \begin{cases}
\ \sigma\_1 & \text{if } \ \sigma\_1 \ge 0 \\
0 & \text{if } \ \sigma\_1 < 0
\end{cases}
\tag{27}
$$

and the maximum compressive stress as:

$$
\sigma\_{\mathfrak{c}} = \begin{cases}
0 & \text{if } \sigma\_3 \ge 0
\end{cases}
\tag{28}
$$

The fluctuations of stresses throughout the cycle can be measured by the variable:

$$
\sigma\_f = \sigma\_1^M - \sigma\_3^m \, \tag{29}
$$

where the superscripts *M* and *m* follow the definition given in Equation (26). The absolute maximum principal stress (AMP *σ*) and strain (AMP*ε*) are defined analogously as:

$$\text{AMP}\sigma = \text{Max}\left\{|\sigma\_1|, |\sigma\_3|\right\} \tag{30a}$$

$$\text{AMP}\varepsilon = \text{Max}\left\{ |\varepsilon\_1|, |\varepsilon\_3| \right\}. \tag{30b}$$

The von Mises ( *σ*vonMises) and Tresca ( *σ*Tresca) stresses (in metallic materials, these variables are used in yielding criteria, which are not applicable to bone, but they can be regarded as well as a measure of the stress intensity) as well as the hydrostatic stress ( *<sup>σ</sup>o*) and volumetric strain (*<sup>ε</sup>o*) are also proposed as predictor variables:

$$
\sigma\_{\text{vonMies}} = \sqrt{\frac{(\sigma\_1 - \sigma\_2)^2 + (\sigma\_1 - \sigma\_3)^2 + (\sigma\_2 - \sigma\_3)^2}{2}} \tag{31a}
$$

$$
\sigma\_{\text{Tresca}} = \frac{\sigma\_1 - \sigma\_3}{2} \tag{31b}
$$

$$
\sigma\_o = \frac{\sigma\_1 + \sigma\_2 + \sigma\_3}{3} \tag{31c}
$$

$$
\varepsilon\_{\mathcal{O}} = \varepsilon\_1 + \varepsilon\_2 + \varepsilon\_3. \tag{31d}
$$

Finally, the SED, which has been extensively used as a mechanical stimulus in bone remodelling algorithms, is also proposed as a predictor. SED is given by the following expression in terms of the stress ( *σ*) and strain (*ε*) tensors or their components ( *<sup>σ</sup>ij*, *<sup>ε</sup>ij*):

$$\text{SED} = \frac{1}{2}\boldsymbol{\sigma} : \boldsymbol{\mathfrak{e}} = \frac{1}{2}\,\sigma\_{ij}\,\varepsilon\_{ij}.\tag{32}$$

Beaupré et al. [3] defined the mechanical stimulus that drives bone remodelling, Ψ, as a combination of two factors: the SED and the number of cycles of each activity. Furthermore, these authors proposed to superimpose the effect of all the activities, *i*, performed by the individual during one day, by weighting the SED of each activity, SED*i*, with the corresponding number of cycles, n*i*. The interested reader can consult the details in [3]. In the following, we will assume that only the most representative activity is carried out daily and that the number of cycles is constant. In that case, there is a linear relationship between Ψ and SED [3]:

$$
\Psi = k \cdot \text{SED.} \tag{33}
$$

In other words, we can identify Ψ and SED for correlation purposes. More importantly, Beaupré et al. proposed to take into account the porosity of the tissue to redefine the mechanical stimulus. Thus, Ψ represents the mechanical stimulus at the continuum (or macroscopic) level. SED can be calculated through FEM and if the constant *k* in Equation (33) is known, the mechanical stimulus at the continuum level, Ψ, can also be evaluated for each element of the mesh. Beaupré et al. hypothesized that this mechanical stimulus must be sensed by the existing tissue within the element in order to produce a remodelling response. Therefore, Ψ can be distributed among the tissue existing in the element through porosity, analogously to what localisation procedures do in multiscale approaches, i.e., allowing to move from the macro to the micro scale. To this end, those authors proposed the following mechanical stimulus at the tissue (or microscopic) level [3]:

$$\Psi\_t = \frac{\Psi}{(1-p)^n} \tag{34}$$

where *p* ∈ [0, 1] is the porosity and *n*=2 is the exponent they used [3]. In this way, if the porosity of one element is close to 1, the mechanical stimulus must be distributed among the little existing tissue and this will be heavily overloaded. Later, we will analyse the effect of the exponent *n*. As stated before, the linear, quadratic and power Pearson coefficients were also evaluated for the SED. The rationale for this is based on the theoretical dependence of SED upon density, which can be deduced from previous works found in the literature [2,35]. In the particular case of a uniaxial stress-state, the stress tensor is *σ* = *σ ei* ⊗*ei*, with *ei* being the loading direction and *σ* the applied stress. In such case, and assuming a linearly elastic and isotropic behaviour for bone, the SED can be calculated through Equation (32) as:

$$\text{SED}^\* = \frac{1}{2}\sigma \,\varepsilon = \frac{1}{2} \, E \,\varepsilon^2,\tag{35}$$

where *ε* is the strain tensor and *ε* is the strain in the loading direction. The asterisk has been added to highlight that this expression corresponds to a particular case. Additionally, a typical power correlation between the Young's modulus and the apparent bone density can be assumed, for example Equations (5a) and (7), which would read:

$$E = B\,\rho^{\mathcal{G}},\tag{36}$$

where *B* and *β* are constants. In this case, Ψ can be rewritten using Equations (33) and (35) as:

$$
\Psi^\* = k \frac{1}{2} B \,\varepsilon^2 \,\rho^\circledast = K \,\rho^\circledast \,\tag{37}
$$

where the constants preceding the factor *ρβ* were grouped in a new constant *K*. The righthand side of Equation (37) is only valid if bone is subjected to a constant strain, which would be in accordance with the Mechanostat Theory. Under all these assumptions, Ψ could be related to the bone apparent density through a power law. Recalling Equation (34) and given that (1 − *p*) is equal to the bone volume fraction, which is proportional to the bone density, the mechanical stimulus at the tissue level, Ψ*t*, could also be related to bone density through a power law. For this reason, we will check if such a power correlation between apparent density and Ψ (or equivalently SED) or Ψ*t* is suitable.
