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

Abstract

Bone remodelling models (BRM) are often used to estimate the density distribution in bones from the loads they are subjected to. BRM define a relationship between a certain variable measuring the mechanical stimulus at each bone site and either the local density or the local variation of density. This agrees with the Mechanostat Theory, which establishes that overloaded bones increase their density, while disused bones tend to decrease their density. Many variables have been proposed as mechanical stimuli, with stress or strain energy density (SED) being some of the most common. Yet, no compelling reason has been given to justify the choice of any of these variables. This work proposes a set of variables derived from the local stress and strain tensors as candidates for mechanical stimuli; then, this work correlates them to the density in the femur of one individual. The stress and strain tensors were obtained from a FE model and the density was obtained from a CT-scan, both belonging to the same individual. The variables that best correlate with density are the stresses. Strains are quite uniform across the femur and very poorly correlated with density, as is the SED, which is, therefore, not a good variable to measure the mechanical stimulus.

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,2,3,4,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:

• Assign an initial uniform density distribution to the bone under study.
• Apply the loads and boundary conditions to a Finite Element (FE) model of the bone and calculate the stresses/strains at every point of the mesh.
• Apply the BRM to relate the stress/strain state with the new density [9] or with the change in density [1,2,3,4,5], and then update the density at each point.
• Update the stiffness tensor at each point based on the new density. Go back to step 2 to start a new iteration.

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, $\dot{\rho }$:

$$\dot{\rho }=f\left(S\right)\mathrm{Evolutionary}\mathrm{model}$$

(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\left(S\right)\mathrm{Non-evolutionary}\mathrm{model}.$$

(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,13,14,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,2,3,4,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 Section 2.4, Section 2.5 and Section 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 Section 2.7 and Section 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).

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 segment 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 force-sharing 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,24,25,26]:

$${\rho }_e=A+B·g_e,$$

(3)

where $g_e$ and ${\rho }_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${}^3$ [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 ${\rho }_e=0$ for $g_e=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 }_e=2.1·\frac{g_e-GLT}{255-GLT}.$$

(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/cm${}^3$), 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 $\nu$ = 0.32 [28].

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]:

$$\begin{array}{c}\hfill E_{\mathrm{Jacobs}}\left(\mathrm{MPa}\right)=\left\{\begin{array}{cc}2014{\rho }^{2.5}\hfill & \mathrm{if} \rho \le 1.2 \mathrm{g}/{\mathrm{cm}}^3\hfill \\ 1763{\rho }^{3.2}\hfill & \mathrm{if} \rho >1.2 \mathrm{g}/{\mathrm{cm}}^3\hfill \end{array}\right.\end{array}$$

(5a)

$$\begin{array}{c}\hfill \nu =\left\{\begin{array}{cc}0.2\hfill & \mathrm{if} \rho \le 1.2 \mathrm{g}/{\mathrm{cm}}^3\hfill \\ 0.32\hfill & \mathrm{if} \rho >1.2 \mathrm{g}/{\mathrm{cm}}^3\hfill \end{array}\right..\end{array}$$

(5b)

Hernandez et al. [29] proposed another correlation for E, based on the ash fraction, $\alpha$ (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=\frac{\rho }{\widehat{\rho }}$, with $\rho$ and $\widehat{\rho }$ being the apparent density and the tissue density, respectively, then the correlation proposed by Hernandez et al. reads:

$$E\left(\mathrm{MPa}\right)=\mathrm{84,370}{\left(\frac{\rho }{\widehat{\rho }}\right)}^{2.58}{\alpha }^{2.74}.$$

(6)

If typical values of ash fraction $\alpha =0.68$ and tissue density $\widehat{\rho }=2.31$ ${\mathrm{g}/\mathrm{cm}}^3$ [29] are used, Equation (6) becomes:

$$E_{\mathrm{Hernandez}}\left(\mathrm{MPa}\right)=3388{\rho }^{2.58},$$

(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 $\nu =0.3$ was used in conjunction with Equation (7) in model IsoH.

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, ${\mathrm{\Psi }}_t$, and the bone formation/resorption rate, $\dot{r}$ (see Figure 4). The reference value of the stimulus, ${\mathrm{\Psi }}_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 $\dot{r}$ determines the temporal evolution of apparent density, $\dot{\rho }$. In turn, apparent density determines the variation of the Young’s modulus through correlations such as (5a) or (7).

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, ${\mathbf{C}}_{\mathbf{0}}$, and damage:

$$\mathbf{C}=(1-D){\mathbf{C}}_{\mathbf{0}},$$

(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 $\mathbf{D}$ and the resulting material is orthotropic, with the principal directions of orthotropy aligning with the principal axes of the damage tensor $\mathbf{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, $\mathbf{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,$$

(9)

with $\mathbf{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, $h_i$:

$$\begin{array}{c}{\displaystyle \frac{1}{E_i}=\frac{1}{\widehat{E}}\frac{1}{h_i^4}}\hfill \\ {\displaystyle -\frac{{\nu }_{ij}}{E_i}=-\frac{\widehat{\nu }}{\widehat{E}}\frac{1}{h_i^2{h}_j^2}}\hfill \\ {\displaystyle \frac{1}{G_{ij}}=\frac{1+\widehat{\nu }}{\widehat{E}}\frac{1}{h_i^2{h}_j^2},}\hfill \end{array}$$

(10)

where $\widehat{E}$ and $\widehat{\nu }$ 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 $\rho$ equal to the density of the bone matrix, $\widehat{\rho }$, which is assumed constant. In this particular modelAnisoH, we assumed $\widehat{\rho }=2.1$ g/cm${}^3$ and chose the Hernandez correlation (7), together with $\nu =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 }{\widehat{\rho }}\right)}^{\beta /4}{A}^{1/4}{\widehat{\mathbf{H}}}^{1/2},$$

(11)

where $\widehat{\mathbf{H}}$ is the normalised fabric tensor. This tensor is normalised by imposing $\det (\widehat{\mathbf{H}})=1$ in order to account only for the directionality of the pores. Additionally, $\beta$ 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 $\rho /\widehat{\rho }$ 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, $\mathbf{Y}$:

$$\mathbf{Y}=2\left\{2\widehat{G}\mathrm{sym}\left[\left(\mathbf{H}\epsilon \mathbf{H}\right)\left(\mathbf{H}\epsilon \right)\right]+\widehat{\lambda }\mathrm{tr}\left({\mathbf{H}}^2\epsilon \right)\mathrm{sym}\left(\mathbf{H}\epsilon \right)\right\},$$

(12)

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

$$\mathbf{J}=\frac{1-\omega }{3}\mathrm{tr}\left(\mathbf{Y}\right)\mathbf{1}+\omega \mathrm{dev}\left(\mathbf{Y}\right),$$

(13)

where $\mathbf{1}$ is the identity tensor, $\mathrm{dev}(•)$ represents the deviatoric part of a tensor and the anisotropy factor, $\omega$, weights the importance of the anisotropy of the stimulus in the model. This factor ranges from $\omega$ = 0, which means that the model only depends on the isotropic component of the stimulus, to $\omega$ = 1, which produces the maximum level of anisotropy. The same value used by Doblaré and García [4] was used here ($\omega$ = 0.1). Two functions $g_r$ and $g_f$ are proposed to establish the remodelling criteria. These functions depend on the stimulus $\mathbf{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, ${\mathrm{\Psi }}_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}{ccc}{g}_f(\mathbf{J},{\mathrm{\Psi }}_t^{*},\mathrm{w})\le 0& g_r(\mathbf{J},{\mathrm{\Psi }}_t^{*},\mathrm{w})>0& \mathrm{resorption};\\ g_r(\mathbf{J},{\mathrm{\Psi }}_t^{*},\mathrm{w})\le 0& g_f(\mathbf{J},{\mathrm{\Psi }}_t^{*},\mathrm{w})>0& \mathrm{formation};\\ g_f(\mathbf{J},{\mathrm{\Psi }}_t^{*},\mathrm{w})\le 0& g_r(\mathbf{J},{\mathrm{\Psi }}_t^{*},\mathrm{w})\le 0& \mathrm{lazy}\mathrm{zone}.\end{array}$$

(14)

Based on the fulfilment of the corresponding criterion, the variation of the fabric tensor $\mathbf{H}$, that accounts for the variation of anisotropy (through tensor $\widehat{\mathbf{H}}$) and porosity (see Equation (11)), is provided by:

$$\begin{array}{cc}{\displaystyle \dot{\mathbf{H}}=k_1\frac{\widehat{\rho }}{\rho }{\mathbf{J}}^{-3}\widehat{\omega }}& \mathrm{resorption};\hfill \\ {\displaystyle \dot{\mathbf{H}}=k_2\frac{\widehat{\rho }}{\rho }\mathbf{J}\widehat{\omega }}& \mathrm{formation};\hfill \\ \dot{\mathbf{H}}=0& \mathrm{lazy}\mathrm{zone};\hfill \end{array}$$

(15)

where the tensor $\widehat{\omega }$ is introduced to simplify the expression, as follows:

$$\widehat{\omega }=\frac{1-2\omega }{3}\mathbf{1}\otimes \mathbf{1}+\omega \mathbf{I}$$

(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 $\dot{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 $\dot{\mathbf{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):

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

(17)

and given that $\widehat{\mathbf{H}}$ must remain normalised ($\det (\widehat{\mathbf{H}})=1$), we finally adopted:

$$\dot{\widehat{\mathbf{H}}}=c{\widehat{\mathbf{H}}}^{1/2}\dot{\mathbf{H}},$$

(18)

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

$$\widehat{\mathbf{H}}\left(t_{j+1}\right)=\widehat{\mathbf{H}}\left(t_j\right)+c{\widehat{\mathbf{H}}}^{1/2}\left(t_j\right)\dot{\mathbf{H}}\left(t_j\right)\Delta t,$$

(19)

where $t_{j+1}$ and $t_j$ are two consecutive integration steps, the time step $\Delta t=1$ is chosen and c is the constant necessary to enforce the condition det $\widehat{\mathbf{H}}\left({\mathbf{t}}_{\mathbf{j}+\mathbf{1}}\right)=1$. The simulation started from an initially isotropic material ($\widehat{\mathbf{H}}\left(t_1\right)$ 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.:

$${\displaystyle \frac{{\sum }_{e=1}^n{\widehat{\mathbf{H}}}_e\left(t_{j+1}\right):{\widehat{\mathbf{H}}}_e\left(t_{j+1}\right)-{\sum }_{e=1}^n{\widehat{\mathbf{H}}}_e\left(t_j\right):{\widehat{\mathbf{H}}}_e\left(t_j\right)}{{\sum }_{e=1}^n{\widehat{\mathbf{H}}}_e\left(t_j\right):{\widehat{\mathbf{H}}}_e\left(t_j\right)}<0.001.}$$

(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.

The peaks of the stimulus coincide with the maximum of the formation criterion function, $g_f$, which is proportional to the stimulus. Those peaks are termed here, $g_f^M$. Since the resorption criterion function, $g_r$, is inversely proportional to the stimulus [4], the local minimum of this function, $g_r^m$, also coincides with the peaks of the stimulus (it must be noted that $g_f^M$ and $g_r^m$ do not necessarily coincide. They are simply reached at the same instant). Analogously, the valleys of the stimulus coincide with the minimum of $g_f$ and the maximum of $g_r$, termed here $g_f^m$ and $g_r^M$, 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 $g_f$ and $g_r$) at every point of the FE mesh throughout the cycle, in order to capture the peaks $g_f^M$ and $g_r^m$ 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}{cc}{g}_f^M>0\hfill & \mathrm{formation};\hfill \\ g_r^m\le 0g_f^M\le 0\hfill & \mathrm{lazy}\mathrm{zone};\hfill \\ g_f^M\le 0g_r^m>0\hfill & \mathrm{resorption}.\hfill \end{array}$$

(21)

It is important to highlight that $g_f^M$ and $g_r^m$ 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 $\mathsf{\Omega }\subset {\mathbb{R}}^3$ be an open bounded domain and $\Gamma =\partial \mathsf{\Omega }$ be its boundary, assumed to be Lipschitz continuous and divided into two disjoint parts ${\Gamma }_D$ and ${\Gamma }_N$ where Dirichlet and Neumann boundary conditions are applied, respectively. We denote, by $\mathbf{x}=x_i{\mathbf{e}}_i$, a generic point of $\mathsf{\Omega }$ and $\mathbf{n}\left(\mathbf{x}\right)=n_i{\mathbf{e}}_i$ as the outward unit normal vector to $\Gamma$ at a point $\mathbf{x}$. The Einstein summation notation was adopted and a Cartesian basis ${\mathbf{e}}_i$ ($i=1,2,3$) can be used without loss of generality.

Let $\mathbf{u}=u_i{\mathbf{e}}_i$, $\mathit{\sigma }={\sigma }_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ and $\mathit{\epsilon }\left(\mathbf{u}\right)={\epsilon }_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ denote the displacement field, the stress tensor and the linearised strain tensor, respectively. ${\mathbf{e}}_i\otimes {\mathbf{e}}_j$ ($i,j=1,2,3$) represents the Cartesian tensorial basis. Let $\mathbf{b}=b_i{\mathbf{e}}_i$ denote the known vector of body forces, $\mathbf{t}=t_i{\mathbf{e}}_i$ the known vector of surface traction forces at ${\Gamma }_N$ and $\mathit{\delta }={\delta }_i{\mathbf{e}}_i$ the known displacements at ${\Gamma }_D$.

The Momentum Conservation Principle states:

$${\sigma }_{ij,j}+b_i=\rho {\ddot{u}}_i,$$

(22)

where $,j$ denotes the partial derivative $\partial /\partial x_j$ and $\ddot{(•)}$ 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:

$${\displaystyle {\epsilon }_{ij}\left(\mathbf{u}\right)=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)i,j=1,2,3}$$

(23)

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

$${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}i,j,k,l=1,2,3,$$

(24)

where $C_{ijkl}$ are the components of the fourth-rank stiffness tensor $\mathbf{C}$. In the case of an isotropic material, this tensor is completely defined by the Young’s modulus, E, and the Poisson’s ratio, $\nu$, 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 ${\Gamma }_D$ and ${\Gamma }_N$:

$$\mathbf{u}=\mathit{\delta }\mathrm{at}{\Gamma }_D$$

(25a)

$${\sigma }_{ij}{n}_j=t_i\mathrm{at}{\Gamma }_N.$$

(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 ($\rho {\ddot{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 $S_e^i$. Then, the following maximum, minimum and amplitude were defined to represent its evolution throughout the gait cycle:

$$S_e^M=\underset{i}{\mathrm{Max}}\left(S_e^i\right)$$

(26a)

$$S_e^m=\underset{i}{\mathrm{Min}}\left(S_e^i\right)$$

(26b)

$$S_e^A=S_e^M-S_e^m.$$

(26c)

Note that the definition of $S_e^M$ 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_e^M$, 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_e^A$.

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${}_j^M$, between the maximum throughout the cycle $S_e^M$ and the apparent density ${\rho }_e$.
• R${}_j^A$, between the amplitude throughout the cycle $S_e^A$ and the apparent density ${\rho }_e$.

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, ${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3$, and analogously, ${\epsilon }_1\ge {\epsilon }_2\ge {\epsilon }_3$. The maximum and minimum principal stresses (${\sigma }_1$, ${\sigma }_3$) and strains (${\epsilon }_1$, ${\epsilon }_3$) are proposed as predictor variables. The maximum tensile stress is defined as:

$${\displaystyle {\sigma }_t=\left\{\begin{array}{cc}{\sigma }_1\hfill & \mathrm{if}{\sigma }_1\ge 0\hfill \\ 0\hfill & \mathrm{if}{\sigma }_1<0\hfill \end{array}\right.}$$

(27)

and the maximum compressive stress as:

$${\displaystyle {\sigma }_c=\left\{\begin{array}{cc}-{\sigma }_3\hfill & \mathrm{if}{\sigma }_3<0\hfill \\ 0\hfill & \mathrm{if}{\sigma }_3\ge 0\hfill \end{array}\right..}$$

(28)

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

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

(29)

where the superscripts M and m follow the definition given in Equation (26). The absolute maximum principal stress (AMP$\sigma$) and strain (AMP$\epsilon$) are defined analogously as:

$${\displaystyle \mathrm{AMP}\sigma =\mathrm{Max}\left\{|{\sigma }_1|,|{\sigma }_3|\right\}}$$

(30a)

$${\displaystyle \mathrm{AMP}\epsilon =\mathrm{Max}\left\{|{\epsilon }_1|,|{\epsilon }_3|\right\}.}$$

(30b)

The von Mises (${\sigma }_{\mathrm{vonMises}}$) and Tresca (${\sigma }_{\mathrm{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 (${\sigma }_o$) and volumetric strain (${\epsilon }_o$) are also proposed as predictor variables:

$${\displaystyle {\sigma }_{\mathrm{vonMises}}=\sqrt{\frac{{({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_1-{\sigma }_3)}^2+{({\sigma }_2-{\sigma }_3)}^2}{2}}}$$

(31a)

$${\displaystyle {\sigma }_{\mathrm{Tresca}}=\frac{{\sigma }_1-{\sigma }_3}{2}}$$

(31b)

$${\displaystyle {\sigma }_o=\frac{{\sigma }_1+{\sigma }_2+{\sigma }_3}{3}}$$

(31c)

$${\displaystyle {\epsilon }_o={\epsilon }_1+{\epsilon }_2+{\epsilon }_3.}$$

(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 ($\mathit{\sigma }$) and strain ($\mathit{\epsilon }$) tensors or their components (${\sigma }_{ij}$, ${\epsilon }_{ij}$):

$${\displaystyle \mathrm{SED}=\frac{1}{2}\mathit{\sigma }:\mathit{\epsilon }=\frac{1}{2}{\sigma }_{ij}{\epsilon }_{ij}.}$$

(32)

Beaupré et al. [3] defined the mechanical stimulus that drives bone remodelling, $\mathrm{\Psi }$, 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 $\mathrm{\Psi }$ and SED [3]:

$$\mathrm{\Psi }=k·\mathrm{SED}.$$

(33)

In other words, we can identify $\mathrm{\Psi }$ 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, $\mathrm{\Psi }$ 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, $\mathrm{\Psi }$, 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, $\mathrm{\Psi }$ 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]:

$${\displaystyle {\mathrm{\Psi }}_t=\frac{\mathrm{\Psi }}{{(1-p)}^n},}$$

(34)

where $p\in [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 $\mathit{\sigma }=\sigma {\mathit{e}}_i\otimes {\mathit{e}}_i$, with ${\mathit{e}}_i$ being the loading direction and $\sigma$ 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:

$${\displaystyle {\mathrm{SED}}^{*}=\frac{1}{2}\sigma \epsilon =\frac{1}{2}E{\epsilon }^2,}$$

(35)

where $\mathit{\epsilon }$ is the strain tensor and $\epsilon$ 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 }^{\beta },$$

(36)

where B and $\beta$ are constants. In this case, $\mathrm{\Psi }$ can be rewritten using Equations (33) and (35) as:

$${\displaystyle {\mathrm{\Psi }}^{*}=k\frac{1}{2}B{\epsilon }^2{\rho }^{\beta }=K{\rho }^{\beta },}$$

(37)

where the constants preceding the factor ${\rho }^{\beta }$ were grouped in a new constant K. The right-hand 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, $\mathrm{\Psi }$ 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, ${\mathrm{\Psi }}_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 $\mathrm{\Psi }$ (or equivalently SED) or ${\mathrm{\Psi }}_t$ is suitable.

3. Results

The weighted Spearman correlation coefficients between bone density and the predictors proposed here are shown in

Table 1. Weighted Spearman correlation coefficients obtained with the constitutive model, AnisoH, using C3D4 elements and for different values of GLT: for the maximum variable throughout the cycle, R${}_S^M$, and for the amplitude throughout the cycle, R${}_S^A$. The predictors analysed are: the maximum principal stress (${\sigma }_1$), the minimum principal stress (${\sigma }_3$), the absolute maximum principal stress (AMP$\sigma$), the corresponding strains (${\epsilon }_1$, ${\epsilon }_3$ and AMP$\epsilon$), the strain energy density (SED or mechanical stimulus at the continuum level $\mathrm{\Psi }$), the mechanical stimulus at the tissue level (${\mathrm{\Psi }}_t$), the maximum tensile and compressive stresses (${\sigma }_t$ and ${\sigma }_c$, respectively), the fluctuation of stresses (${\sigma }_f$), the von Mises and Tresca stresses, the hydrostatic stress (${\sigma }_o$) and the volumetric strain (${\epsilon }_o$). Values higher than 0.9 are highlighted in boldface; negative values are in red.

Predictor	GLT = 70		GLT = 80		GLT = 90
	R${}_{\mathit{S}}^{\mathit{M}}$	R${}_{\mathit{S}}^{\mathit{A}}$	R${}_{\mathit{S}}^{\mathit{M}}$	R${}_{\mathit{S}}^{\mathit{A}}$	R${}_{\mathit{S}}^{\mathit{M}}$	R${}_{\mathit{S}}^{\mathit{A}}$
${\sigma }_1$	0.871	0.880	0.905	0.912	0.936	0.942
${\sigma }_3$	−0.040	0.872	−0.047	0.897	−0.045	0.926
${\sigma }_t$	0.871	0.871	0.905	0.905	0.936	0.936
${\sigma }_c$	0.868	0.868	0.894	0.894	0.923	0.923
${\sigma }_f$	0.904	0.904	0.927	0.927	0.951	0.951
AMP$\sigma$	0.896	0.896	0.921	0.920	0.945	0.944
${\sigma }_{\mathrm{vonMises}}$	0.898	0.898	0.922	0.922	0.948	0.948
${\sigma }_{\mathrm{Tresca}}$	0.898	0.898	0.922	0.922	0.948	0.947
${\sigma }_o$	0.701	0.886	0.739	0.911	0.782	0.938
${\epsilon }_1$	−0.344	−0.340	−0.431	−0.429	−0.503	−0.502
${\epsilon }_3$	0.410	−0.355	0.432	−0.445	0.457	−0.504
AMP$\epsilon$	−0.328	−0.325	−0.427	−0.425	−0.499	−0.498
$\mathrm{\Psi }$	0.698	0.698	0.721	0.721	0.772	0.772
${\mathrm{\Psi }}_t$ for n = 2	−0.048	−0.048	−0.092	−0.092	−0.097	−0.099
${\epsilon }_o$	−0.202	−0.303	−0.255	−0.468	−0.327	−0.571

for the constitutive model AnisoH. As stated before, three different values were used for the grey level threshold (GLT) used in Equation (4), thus, leading to three different FE models (see Figure 2). The correlation coefficients are given in the three cases for comparison.

The fact that the Spearman correlation coefficient is non-parametric makes it more appropriate to evaluate the correlation between density and the variables, as it implies no assumption on the type of relationship. It simply establishes if there is a concordance between those points having the highest density and those having the highest values of a certain predictor.

It can be seen that most of the stress magnitudes are highly correlated with the density except for the peak of the minimum principal stress, ${\sigma }_3^M$, for obvious reasons, as the sign of ${\sigma }_3$ (usually negative) is considered in the calculus of this peak. Therefore, ${\sigma }_3^M$ usually corresponds to the lowest absolute value throughout the cycle. The amplitude, ${\sigma }_3^A$, is better correlated with the density as it usually measures the range of the compressive stress.

The strain magnitudes are poorly correlated with the density, in some cases with negative coefficients. SED (or $\mathrm{\Psi }$) is only moderately correlated with the density and not correlated at all if it is corrected to account for the porosity (${\mathrm{\Psi }}_t$).

Regarding the influence of GLT, it can be seen that though the values of R are different, the same trend is observed in the three cases. In fact, if we ordered the predictors based on R, the same order would result for the three GLT.

The weighted Pearson correlation coefficients are shown in

Table 2. Weighted Pearson correlation coefficients obtained with the constitutive model, AnisoH, using C3D4 elements and for different values of GLT: for the maximum variable throughout the cycle, R${}_P^M$, and for the amplitude throughout the cycle, R${}_P^A$. The predictors analysed in this case are: the absolute maximum principal stress (AMP$\sigma$), the maximum tensile and compressive stresses (${\sigma }_t$ and ${\sigma }_c$, respectively), the fluctuation of stresses (${\sigma }_f$) and the mechanical stimulus at the continuum level ($\mathrm{\Psi }$) and at the tissue level (${\mathrm{\Psi }}_t$) for two values of n. Values higher than 0.8 are highlighted in boldface; negative values are in red.

Predictor	Type	GLT = 70		GLT = 80		GLT = 90
		R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$	R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$	R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$
AMP$\sigma$	Linear	0.773	0.772	0.793	0.791	0.808	0.806
${\sigma }_t$	Linear	0.746	0.745	0.761	0.760	0.770	0.769
${\sigma }_c$	Linear	0.671	0.672	0.684	0.684	0.695	0.695
${\sigma }_f$	Linear	0.810	0.810	0.825	0.825	0.838	0.838
$\mathrm{\Psi }$	Linear	0.012	0.012	0.012	0.012	0.012	0.012
${\mathrm{\Psi }}_t$ for n = 1	Linear	−0.016	−0.016	−0.016	−0.016	−0.013	−0.013
${\mathrm{\Psi }}_t$ for n = 2	Linear	−0.023	−0.023	−0.026	−0.026	−0.023	−0.023
AMP$\sigma$	Quadratic	0.784	0.782	0.804	0.803	0.818	0.817
${\sigma }_t$	Quadratic	0.762	0.761	0.778	0.778	0.787	0.786
${\sigma }_c$	Quadratic	0.678	0.679	0.690	0.691	0.700	0.700
${\sigma }_f$	Quadratic	0.823	0.823	0.839	0.839	0.850	0.850
$\mathrm{\Psi }$	Quadratic	0.026	0.026	0.026	0.026	0.029	0.029
${\mathrm{\Psi }}_t$ for n = 1	Quadratic	0.017	0.017	0.017	0.017	0.018	0.018
${\mathrm{\Psi }}_t$ for n = 2	Quadratic	0.021	0.021	0.026	0.026	0.023	0.023
AMP$\sigma$	Power	0.780	0.778	0.799	0.798	0.813	0.811
${\sigma }_t$	Power	0.761	0.760	0.777	0.776	0.785	0.784
${\sigma }_c$	Power	0.674	0.675	0.683	0.684	0.690	0.691
${\sigma }_f$	Power	0.820	0.820	0.835	0.835	0.845	0.845
$\mathrm{\Psi }$	Power	0.006	0.006	0.005	0.005	0.005	0.005
${\mathrm{\Psi }}_t$ for n = 1	Power	−0.016	−0.016	−0.017	−0.017	−0.015	−0.015
${\mathrm{\Psi }}_t$ for n = 2	Power	−0.019	−0.019	−0.022	−0.020	−0.015	−0.015

for the constitutive model, AnisoH. Only some of the variables previously analysed in Table 1—those that are considered more interesting—are studied here; in particular, some of the stress variables that had a higher Spearman coefficient together with the SED at the continuum and at the tissue level for n = 1 and n = 2 (see Equation (34)). The Pearson coefficients are parametric and presuppose a certain relationship between the variables being correlated. Thus, we have tried linear, quadratic and power functions (see Appendix A for details).

Compared to the Spearman, the Pearson correlations have worsened notably as we are forcing them to fit a certain function which is probably not the most appropriate to relate the density with the predictor. Among the three types of functions tested, the quadratic is slightly better, followed by the power and the linear function. The low correlation between SED and density stands out—something that does not improve in the case of SED at the tissue level—for which even negative correlation coefficients were obtained, as in the case of the Spearman coefficients.

Table 3. Weighted Spearman correlation coefficients obtained for GLT = 80, constitutive models AnisoH, IsoH and IsoJ and using C3D4 elements: for the maximum variable throughout the cycle, R${}_S^M$, and for the amplitude throughout the cycle, R${}_S^A$. The predictors analysed are: the maximum principal stress (${\sigma }_1$), the minimum principal stress (${\sigma }_3$), the absolute maximum principal stress (AMP$\sigma$), the corresponding strains (${\epsilon }_1$, ${\epsilon }_3$ and AMP$\epsilon$), the strain energy density (SED or mechanical stimulus at the continuum level $\mathrm{\Psi }$), the mechanical stimulus at the tissue level (${\mathrm{\Psi }}_t$), the maximum tensile and compressive stresses (${\sigma }_t$ and ${\sigma }_c$, respectively), the fluctuation of stresses (${\sigma }_f$), the von Mises and Tresca stresses, the hydrostatic stress (${\sigma }_o$) and the volumetric strain (${\epsilon }_o$). Values higher than 0.9 are highlighted in boldface; negative values are in red.

Predictor	AnisoH		IsoH		IsoJ
	R${}_{\mathit{S}}^{\mathit{M}}$	R${}_{\mathit{S}}^{\mathit{A}}$	R${}_{\mathit{S}}^{\mathit{M}}$	R${}_{\mathit{S}}^{\mathit{A}}$	R${}_{\mathit{S}}^{\mathit{M}}$	R${}_{\mathit{S}}^{\mathit{A}}$
${\sigma }_1$	0.905	0.912	0.906	0.913	0.905	0.913
${\sigma }_3$	−0.047	0.897	−0.058	0.899	−0.059	0.895
${\sigma }_t$	0.905	0.905	0.906	0.906	0.905	0.905
${\sigma }_c$	0.894	0.894	0.895	0.895	0.892	0.892
${\sigma }_f$	0.927	0.927	0.928	0.928	0.925	0.925
AMP$\sigma$	0.921	0.920	0.922	0.921	0.920	0.920
${\sigma }_{\mathrm{vonMises}}$	0.922	0.922	0.924	0.924	0.921	0.920
${\sigma }_{\mathrm{Tresca}}$	0.922	0.922	0.923	0.923	0.920	0.920
${\sigma }_o$	0.739	0.911	0.742	0.913	0.739	0.911
${\epsilon }_1$	−0.431	−0.429	−0.431	−0.429	−0.443	−0.441
${\epsilon }_3$	0.432	−0.445	0.431	−0.443	0.429	−0.458
AMP$\epsilon$	−0.427	−0.425	−0.425	−0.423	−0.440	−0.438
$\mathrm{\Psi }$	0.721	0.721	0.722	0.721	0.708	0.708
${\mathrm{\Psi }}_t$ for n = 2	−0.092	−0.092	−0.092	−0.092	−0.114	−0.114
${\epsilon }_o$	−0.255	−0.468	−0.254	−0.464	−0.247	−0.489

and

Table 4. Weighted Pearson correlation coefficients obtained for GLT = 80, constitutive models AnisoH, IsoH and IsoJ and using C3D4 elements: for the maximum variable throughout the cycle, R${}_P^M$, and for the amplitude throughout the cycle, R${}_P^A$. The predictors analysed in this case are: the absolute maximum principal stress (AMP$\sigma$), the maximum tensile and compressive stresses (${\sigma }_t$ and ${\sigma }_c$, respectively), the fluctuation of stresses (${\sigma }_f$), the mechanical stimulus at the continuum level ($\mathrm{\Psi }$) and at the tissue level (${\mathrm{\Psi }}_t$) for two values of n. Values higher than 0.8 are highlighted in boldface; negative values are in red.

Predictor	Type of Correlation	AnisoH		IsoH		IsoJ
		R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$	R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$	R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$
AMP$\sigma$	Linear	0.793	0.791	0.802	0.800	0.814	0.813
${\sigma }_t$	Linear	0.761	0.760	0.767	0.766	0.771	0.770
${\sigma }_c$	Linear	0.684	0.684	0.692	0.692	0.707	0.708
${\sigma }_f$	Linear	0.825	0.825	0.832	0.832	0.840	0.840
$\mathrm{\Psi }$	Linear	0.012	0.012	0.012	0.012	0.015	0.015
${\mathrm{\Psi }}_t$ for n = 1	Linear	−0.016	−0.016	−0.016	−0.016	−0.017	−0.017
${\mathrm{\Psi }}_t$ for n = 2	Linear	−0.026	−0.026	−0.026	−0.026	−0.029	−0.029
AMP$\sigma$	Quadratic	0.804	0.803	0.814	0.812	0.834	0.833
${\sigma }_t$	Quadratic	0.778	0.778	0.785	0.784	0.797	0.796
${\sigma }_c$	Quadratic	0.690	0.691	0.698	0.699	0.719	0.720
${\sigma }_f$	Quadratic	0.839	0.839	0.846	0.846	0.861	0.861
$\mathrm{\Psi }$	Quadratic	0.026	0.026	0.026	0.026	0.030	0.030
${\mathrm{\Psi }}_t$ for n = 1	Quadratic	0.017	0.017	0.017	0.017	0.020	0.020
${\mathrm{\Psi }}_t$ for n = 2	Quadratic	0.026	0.026	0.026	0.026	0.029	0.029
AMP$\sigma$	Power	0.799	0.798	0.809	0.808	0.833	0.832
${\sigma }_t$	Power	0.777	0.776	0.783	0.783	0.797	0.796
${\sigma }_c$	Power	0.683	0.684	0.692	0.693	0.717	0.718
${\sigma }_f$	Power	0.835	0.835	0.842	0.842	0.861	0.861
$\mathrm{\Psi }$	Power	0.005	0.005	0.006	0.006	0.008	0.008
${\mathrm{\Psi }}_t$ for n = 1	Power	−0.017	−0.017	−0.017	−0.017	−0.018	−0.018
${\mathrm{\Psi }}_t$ for n = 2	Power	−0.022	−0.022	−0.022	−0.022	−0.023	−0.023

compare, respectively, the Spearman and Pearson coefficients obtained using the three constitutive models analysed in this work: AnisoH, IsoH and IsoJ. As indicated above, the effect of GLT was not important, and hence, only one case (GLT = 80) was studied. It can be noted that the effect of the constitutive model is negligible on the Spearman correlations and small on the Pearson correlations, at least for the three cases tested here. The biggest difference is obtained for the power fit between the IsoH and IsoJ models, which, in turn, follow a different power correlation between the density and the Young’s modulus, but is not greater than 0.03.

We have also analysed the spatial distribution of the correlations, in particular of the Pearson coefficients (power fit), by assessing separately the correlations for the elements of the proximal, distal and diaphyseal thirds (see

Table 5. Weighted Pearson correlation coefficients (power fit) obtained for the maximum variable throughout the cycle, R${}_S^M$, for GLT = 80, constitutive model IsoH and using C3D4 elements. The predictors analysed are: the fluctuation of stresses (${\sigma }_f$) and the absolute maximum principal stress (AMP$\sigma$). Values higher than 0.8 are highlighted in boldface.

	Proximal	Diaphysis	Distal	Global
${\sigma }_f$	0.718	0.879	0.660	0.844
AMP$\sigma$	0.694	0.862	0.640	0.813

). The aim of this comparison was to investigate if there are regions of the femur where the density is better to the predictors. Given the limited influence of GLT and the constitutive model, we only show the case GLT = 80 and IsoH. Besides, we only compare some of the variables that show a higher correlation (${\sigma }_f$, AMP$\sigma$) and only the coefficients for the maximum variable throughout the cycle, R${}_S^M$. The other stress variables follow the same trend, as well as R${}_S^A$. The comparison of the strain variables is meaningless since they are not correlated with density, as shown previously. It can be noted that the correlation coefficients are high in the diaphysis, significantly worse in the proximal and especially worse in the distal third, influenced by the simplified way the joint reaction forces were modelled. They were applied as concentrated nodal forces, as explained in Section 2.7, rather than as a load distributed over the articular surface, as it actually occurs. This simplification affects the stresses near the articular region and, therefore, the correlations. The hip joint force can be more plausibly applied as a concentrated nodal force since the pressure on that joint spans a narrower region than that on the knee joint. Probably, this makes the correlations be slightly better in the proximal third than in the distal one.

The influence of the mesh (C3D4 vs. C3D10) was analysed by comparing the correlation coefficients in

Table 6. Influence of the FE mesh (C3D4 vs. C3D10) on the weighted Spearman and Pearson (power) correlation coefficients obtained for GLT = 80 and constitutive model IsoH: for the maximum variable throughout the cycle, R${}_P^M$, and for the amplitude throughout the cycle, R${}_P^A$. The same predictors analysed for Spearman and Pearson coefficients are compared here. Values higher than 0.9 (in Spearman) or 0.8 (in Pearson) are highlighted in boldface; negative values are in red.

Predictor	Type of Correlation	C3D4		C3D10
		R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$	R${}_{\mathit{P}}^{\mathit{M}}$	R${}_{\mathit{P}}^{\mathit{A}}$
${\sigma }_1$	Spearman	0.906	0.913	0.904	0.912
${\sigma }_3$	Spearman	−0.058	0.899	−0.236	0.896
${\sigma }_t$	Spearman	0.906	0.906	0.904	0.904
${\sigma }_c$	Spearman	0.895	0.895	0.893	0.894
${\sigma }_f$	Spearman	0.928	0.928	0.930	0.930
AMP$\sigma$	Spearman	0.922	0.921	0.923	0.922
${\sigma }_{\mathrm{vonMises}}$	Spearman	0.924	0.924	0.925	0.925
${\sigma }_{\mathrm{Tresca}}$	Spearman	0.923	0.923	0.925	0.924
${\sigma }_o$	Spearman	0.742	0.913	0.734	0.911
${\epsilon }_1$	Spearman	−0.431	−0.429	−0.499	−0.496
${\epsilon }_3$	Spearman	0.431	−0.443	0.484	−0.495
AMP$\epsilon$	Spearman	−0.425	−0.423	−0.482	−0.480
$\mathrm{\Psi }$	Spearman	0.722	0.721	0.709	0.709
${\mathrm{\Psi }}_t$ for n = 2	Spearman	−0.092	−0.092	−0.182	−0.182
${\epsilon }_o$	Spearman	−0.254	−0.464	−0.287	−0.499
AMP$\sigma$	Pearson (Power)	0.809	0.808	0.874	0.874
${\sigma }_t$	Pearson (Power)	0.783	0.783	0.818	0.818
${\sigma }_c$	Pearson (Power)	0.692	0.693	0.741	0.742
${\sigma }_f$	Pearson (Power)	0.842	0.842	0.902	0.902
$\mathrm{\Psi }$	Pearson (Power)	0.006	0.006	0.011	0.011
${\mathrm{\Psi }}_t$ for n = 1	Pearson (Power)	−0.017	−0.017	−0.019	−0.019
${\mathrm{\Psi }}_t$ for n = 2	Pearson (Power)	−0.022	−0.022	−0.025	−0.025

; in particular, the Spearman and the power Pearson coefficients in the case GLT = 80 and using the constitutive model IsoH. The Pearson coefficients improved moderately with the use of quadratic elements (C3D10), but only for the good predictors, i.e., those variables that are highly correlated with density. The rest of variables, such as $\mathrm{\Psi }$ and the strain magnitudes (not shown), did not improve their correlations. The other types of fit (linear and quadratic) also improved with C3D10, although to a lesser extent. It is noteworthy that the Spearman coefficients were almost identical in both meshes.

Histograms

There are so many points involved in the correlations that the plots of density against the different variables are very difficult to distinguish. Instead, histograms are used to show the percentage of volume occupied by those elements whose values of AMP$\sigma$ or AMP$\epsilon$ are within a given range. The elements of cortical ($\rho >1.2$ g/cm${}^3$) and trabecular bone ($\rho \le 1.2$ g/cm${}^3$) have been separated into two different histograms and the results for the three constitutive models were plotted jointly (see Figure 6). Thus, for example, in the model AnisoH, the elements of trabecular bone whose AMP$\epsilon$ is in the range [200,600] $\mathsf{\mu }\epsilon$ occupy 30% of the total volume of trabecular bone ($\rho \le 1.2$ g/cm${}^3$).

It can be seen that while the strain range is similar for cortical and trabecular bone, the stresses are completely different, with trabecular bone having stresses several orders of magnitude lower than cortical bone. In general, strains are found in a very narrow range, especially in the cortical bone, for which 87% of the volume has AMP$\epsilon$ in the range of 200–600 $\mathsf{\mu }\epsilon$. This is not so evident in trabecular bone, though 51% is still within the range of 200–1000 $\mathsf{\mu }\epsilon$ and 65% is in the range of 200–1400 $\mathsf{\mu }\epsilon$. This is still a narrow range, as overload strains are up to 4000 $\mathsf{\mu }\epsilon$ [36,37].

4. Discussion

A linear function (Equation (3)) was used here to estimate the bone apparent density from the grey level of the CT scan [23,24,25,26]. As in these previous works, the CT scan could not be calibrated, and thus, a threshold was set by identifying the grey level of the diaphyseal canal with bone marrow, that is, with null bone density. Notwithstanding, the greyscale value inside the canal was found to vary in the range of 70–90 and, thus, three different grey level thresholds (GLT) were used and their effect was investigated. The same trend was observed in R for the three values of GLT. Therefore, we can state that the choice of GLT had no effect on the correlation coefficients (see Table 1 and Table 2), and hence, on the overall conclusions of this study. For that reason, the rest of results were compared only for the intermediate GLT = 80.

In general, the variables derived from the strain tensor has a low correlation with density and this can be due to the fact that the strains are concentrated in a narrow range, as deduced from the histograms. On the contrary, the stress-related variables are distributed over a much wider range and they are very closely related to density as the correlation coefficients showed. The high correlation found between density and the von Mises, Tresca, absolute maximum principal stress (AMP$\sigma$) and the fluctuation of stresses (${\sigma }_f$) are noteworthy.

On the other hand, SED is only modestly correlated with density, probably because it depends on the strains, which are very poorly correlated with $\rho$. Since SED also depends on the stress level, the histograms of SED (not shown) are as spread as the histograms of stress; however, in view of the correlations, its variation does not seem to be as coupled to density as the mere variation of stress is. However, if SED seems a modest predictor of density, the mechanical stimulus at the tissue level proposed by Beaupré (${\mathrm{\Psi }}_t$) is even worse, as it yields negative correlation coefficients.

The Spearman coefficient is probably the most simple indicator that two variables are correlated by a monotonic function, since it does not assume any specific relationship as does the Pearson coefficient. If R${}_S$ is positive and high for a certain variable, this means that density increases with that variable, and that a large part of that increase can be explained by the variable, without assessing what specific relationship exists. Hence, that variable can serve as a good predictor of density in a bone remodelling model, though the particular relationship between the density and the specific variable should be further investigated. In this regard, the linear, quadratic or power functions tried for the Pearson coefficients worked only moderately well. Only ${\sigma }_f$ and AMP$\sigma$ showed a value of R${}_P$ slightly over 0.8 (up to 0.861 for ${\sigma }_f$ in the IsoJ model). This means that R${}^2>\sim$ 0.64, i.e., around 64% of the variation of the density, can be explained by ${\sigma }_f$ (up to 74% in the IsoJ model).

Only some of the variables appearing in Table 1 were chosen for assessing its Pearson coefficient, those having a high Spearman coefficient plus SED and its related variables (${\mathrm{\Psi }}_t$), for its common use as a mechanical stimulus in many models of bone remodelling. Obviously, if the Spearman coefficients of SED and its related variables were low, the respective Pearson coefficients could not improve them, but it is worth noting the very low values of R${}_P$ obtained for SED compared to R${}_S$. This would mean that Equation (37) was not very appropriate, probably because of the many assumptions involved in it, viz: uniaxial stress-state, isotropic material and power correlation between Young’s modulus and density.

The distribution of two variables, one representing the stress-state (AMP$\sigma$) and one representing the strain-state (AMP$\epsilon$), was analysed by means of histograms (see Figure 6) that distinguish between cortical ($\rho >1.2$ g/cm${}^3$) and trabecular bone ($\rho \le 1.2$ g/cm${}^3$). These histograms showed that the strains are concentrated in a relatively narrow range in the case of trabecular bone. Around 65% of the volume of trabecular bone was found to be within the range AMP$\epsilon \in$ [200–1400] $\mathsf{\mu }\epsilon$ in the AnisoH model, 62% in IsoH and 53% in the IsoJ. The range was particularly narrow in the case of cortical bone, as about 85% of its total volume was found in the range AMP$\epsilon \in$ [200–600] $\mathsf{\mu }\epsilon$ for all constitutive models.

It is noteworthy that the range of strains obtained here was low compared with the normal strains indicated by the Mechanostat Theory (between 800 and 1200 $\mathsf{\mu }\epsilon$ [10]). Nonetheless, other authors have established a different range of normal strains in the so-called “adapted-window” of the Mechanostat, between 200 and 1500 $\mathsf{\mu }\epsilon$ [36,37]. Yet, the strains seem relatively low for models, AnisoH and IsoH, which are likely overestimating the bone stiffness. The correlation, (5a), used by model IsoJ, is probably more adequate in view of the strains it produces.

In contrast to strains, the stresses are more uniformly distributed and over a much wider range, with the stresses of trabecular bone being one or two orders of magnitude lower than those of cortical bone. This strong relationship between bone density and stress is confirmed by the high Spearman correlation coefficients of most stress variables (see Table 1).

The distributions and correlations of strain and stress variables would confirm that bone is adapted to withstand a constant strain, or at least a strain within a narrow range, in accordance with the Mechanostat Theory, while the local stress-state seems to determine the bone density. This would suggest the use of a strain variable as the mechanical stimulus S in evolutionary BRM (see Equation (1)), such that bone density changes if the strain is out of the normal or adapted range. On the contrary, a stress variable would be more appropriate for the stimulus in a non-evolutionary BRM (see Equation (2)), such that the stress would determine the bone density at a given location. In the case of cyclic loads, such as the one applied here, a variable measuring the fluctuations of stresses throughout the cycle seems the more appropriate mechanical stimulus among those tested here, though the specific relationship between density and stress is yet to be determined. In no case does the SED seem a suitable variable to be used as the mechanical stimulus in BRM.

We have compared three constitutive models (AnisoH, IsoH and IsoJ) in order to check whether the correlation between predictors and density, $S-\rho$, is being forced by the relationship between stiffness and density, $E=F\left(\rho \right)$, in particular for those predictors derived from the stress tensor (for the sake of generality, we should write $\mathbf{C}=F\left(\rho \right)$, with $\mathbf{C}$ being the stiffness tensor, in order to account for anisotropic materials. However, the rationale provided is independent of this distinction and we will continue with $E=F\left(\rho \right)$). It is well known from the literature that F is a monotonically increasing function and given that the stiffer elements tend to withstand higher stress levels, we should expect a function such as $\sigma =H\left(E\right)$ to be monotonically increasing as well. Thus, if we write:

$$\sigma =H\left(F\left(\rho \right)\right),$$

(38)

we can expect a monotonically increasing function relating $\sigma$ and $\rho$, in other words, a positive correlation $S-\rho$ for the predictors derived from the stress tensor. Therefore, it could seem that the function $E=F\left(\rho \right)$ is forcing $S-\rho$ and it is so to some extent, especially for the Spearman correlation coefficients. However, if $E=F\left(\rho \right)$ were the only determinants of the correlation $S-\rho$, this correlation should depend on the constitutive model and it does not—significantly, at least. Indeed, the increase in stiffness from IsoJ to IsoH is quite remarkable (up to 70%, see Figure 3, left) and yet the Pearson correlations do not change much from one model to the other (see Table 4). Moreover, since the ratio of stiffness is non-uniform (see Figure 3, right), a redistribution of stresses could be expected that would cause the stress patterns of the two models not to coincide. The same could be said of the comparison between models AnisoH and IsoH. Such redistribution of stresses would, therefore, affect the Spearman correlation coefficients or the AMP$\sigma$ histograms, but both are almost identical for the three models (see Table 3 and Figure 6). Hence, we can conclude that the assumed constitutive model has no significant effect on the correlations.

However, there is another argument to support that $S-\rho$ is not completely forced by the constitutive behaviour, but only to some extent. Certainly, Equation (38) is incomplete as the stress-state also depends on the loads and boundary conditions:

$$\sigma =H\left(F\left(\rho \right),\mathbf{F}\right),$$

(39)

where $\mathbf{F}$ stands for the set of loads applied to the domain, including the body forces, the surface tractions applied on the boundaries (i.e., Neumann boundary conditions) and the reaction forces (and consequently the Dirichlet boundary conditions). Apart from that, the predictors, S, are obtained as a function of the stress tensor, that we can denote as $S=G\left(\sigma \right)$, thus leading to:

$$S=G\left(H\left(F\left(\rho \right),\mathbf{F}\right)\right).$$

(40)

In view of Equation (40) the predictors, S, need not be predetermined only by the function $F\left(\rho \right)$, not even the Spearman correlation coefficients, as the function G is not necessarily monotonic. This is confirmed by the different Spearman coefficients obtained for different predictors (see Table 1 and Table 3). It should be noted that the same concepts that are behind the G functions can be applied to the predictors derived from the strain tensor or the SED-based predictors.

An example can serve to illustrate the key role played by the loads on the correlations $S-\rho$. We have obtained these correlations for every instant of the gait cycle separately. The stress and strain tensors obtained at each instant i were used to assess $S_e^i$ (recall Section 2.7) and these variables were correlated with ${\rho }_e$ to give R${}_P^i$, a weighted Pearson correlation coefficient for each instant i. The worst coefficient between AMP$\sigma$ and density throughout the cycle was R${}_P^i$ = 0.48 (for the quadratic fit and in the case GLT = 80, constitutive model IsoH and C3D4 elements), which is far from R${}_P^M$ = 0.804. The former is a very poor correlation. (R${}_P^i)$${}^2$ = 0.23, i.e., only about 25% of the variance of density can be explained by AMP${\sigma }^i$ and this is probably because the loads at that instant are not representative of the gait cycle. In fact, the loads at a single instant cannot be representative of the entire cycle on a general basis.

In summary, the function $E=F\left(\rho \right)$ is contributing to obtain the correlations $S-\rho$, at least for stress-based predictors, but only to some extent. There are other factors, not influenced by that function, that play an important role in the correlations; an accurate estimation of the loads is crucial, taking into account the variation of stresses throughout the cycle is also important and an appropriate choice of S (funtion G) is key to predicting density.

The type of element (C3D4 vs. C3D10) had a relative importance on the results, as the linear elements (C3D4) are stiffer than the quadratic ones (C3D10). The latter yielded moderately higher stresses and strains, which were slightly better correlated with density through a power law than those obtained with C3D4 elements (see Table 6). The linear and quadratic fits (not shown) also improved with C3D10, but to a lesser extent. Notwithstanding, the C3D4 mesh had a sufficient element density and was accurate enough as evidenced by the fact that the stress and strain patterns are almost identical in both meshes. Consequently, the Spearman correlation coefficients derived from both are also almost identical. For that reason, the use of C3D4 elements was justified, at least for the mesh density employed in our model.

Among the limitations of this study we must note that only one activity (walking) has been considered, among the many activities that a person can carry out in a normal day. Other activities can load the femur in a different way than walking, thus, affecting the local bone density. This could partially explain the variance of density not explained by our correlations. However, walking is by far the most common activity affecting the lower limb and the most common activity performed by the subject under study. Hence, little influence of other activities can be expected.

The joint loads were applied in a simplified way, as nodal-concentrated loads instead of loads distributed over the articular surface. This approximation is especially significant in the knee reaction force and this could explain the lower correlation coefficients found in the distal third of the femur. Another limitation is that the CT scan could not be calibrated, but we have shown that this fact did not influence the conclusions drawn. Finally, it must be noted that only one individual, a male healthy subject, was studied. As future work, it would be interesting to repeat the study in other cases, including different bones, age, gender, health status, bone pathologies, etc., in order to confirm whether these variables or others alter the dependence of bone density on the proposed predictors.

5. Conclusions

In view of the results of this study we can conclude the following:

• Stresses are highly variable in the femur and very different in cortical and trabecular sites.
• Stresses are strongly correlated with density and seem to determine it. Therefore, certain variables derived from the stress tensor could be good candidates for the mechanical stimulus in non-evolutionary BRM that can be used in DPTBR problems.
• In the case of cyclic loads, the fluctuation of stresses throughout the cycle is the best predictor of density among those tested in this work.
• In contrast to the above, strains are relatively uniform across the femur and remain within the “adapted-window” established by the Mechanostat Theory.
• Strains are very poorly correlated with density and all indications are that bone tends to keep them approximately constant. Therefore, variables derived from the strain tensor could be good candidates as mechanical stimulus in evolutionary BRM, in which density would be forced to change in order to maintain strains within the “adapted-window”.
• Strain energy density does not seem to be a good variable to measure the mechanical stimulus in any case, and even less so if this variable is evaluated at the tissue level by using the porosity correction factor.