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

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/). mathematics Article 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 José Luis Calvo-Gallego 1,* , Fernando Gutiérrez-Millán 1, Joaquín Ojeda 1, María Ángeles Pérez 2 and Javier Martínez-Reina 1,* 1Departmento de Ingeniería Mecánica y Fabricación, Universidad de Sevilla, 41092 Seville, Spain 2 Multiscale in Mechanical and Biological Engineering, Instituto de Investigación en Ingeniería de Aragón (I3A), Instituto de Investigación Sanitaria Aragón (IIS Aragón), University of Zaragoza, 50018 Zaragoza, Spain *Correspondence: [email protected] (J.L.C.-G.); [email protected] (J.M.-R.) 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. Keywords: bone remodelling; mechanical stimulus; correlation; bone density distribution; strain energy density; absolute maximum principal stress; fluctuation of stresses MSC: 92-10 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: 1. Assign an initial uniform density distribution to the bone under study. 2. 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. Mathematics 2022,10, 3367. https://doi.org/10.3390/math10183367 https://www.mdpi.com/journal/mathematics Mathematics 2022,10, 3367 2 of 29 3. Apply the BRM to relate the stress/strain state with the new density [ 9 ] or with the change in density [1–5], and then update the density at each point. 4. 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, ˙ ρ: ˙ ρ=f(S)Evolutionary 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. ρ=g(S)Non-evolutionary 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 – 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 Mathematics 2022,10, 3367 3 of 29 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 4and 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 1and 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). Mathematics 2022,10, 3367 4 of 29 0 1 2 3 4 0 0.5 1 1.5 2 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 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 Material. 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. Mathematics 2022,10, 3367 5 of 29 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]: ρe=A+B·ge, (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/cm3[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: ρe=2.1 ·ge−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 2shows 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 2that 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]. Mathematics 2022,10, 3367 6 of 29 DENSITY (g/cm3) +1.000.10-3 +1.759.10-1 +3.508.10-1 +5.258.10-1 +7.007.10-1 +8.756.10-1 +1.051 +1.225 +1.400 +1.575 +1.750 +1.925 +2.100 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]: EJacobs (MPa) = 2014 ρ2.5 if ρ≤1.2 g/cm3 1763 ρ3.2 if ρ>1.2 g/cm3(5a) ν=0.2 if ρ≤1.2 g/cm3 0.32 if ρ>1.2 g/cm3. (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(MPa) = 84370 ρ ˆ ρ2.58 α2.74. (6) If typical values of ash fraction α= 0.68 and tissue density ˆ ρ= 2.31 g/cm3 [ 29 ] are used, Equation (6) becomes: EHernandez(MPa) = 3388 ρ2.58, (7) Mathematics 2022,10, 3367 7 of 29 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. 0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 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 Mathematics 2022,10, 3367 8 of 29 linear relationship between the mechanical stimulus, Ψt , and the bone formation/resorption rate, ˙ r (see Figure 4). The reference value of the stimulus, Ψ∗ 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). Ψt−Ψ∗ t ˙r 1 c2 1 c1 w w 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, Ψ∗ 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: C= (1−D)C0, (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: D=1−H2, (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: Mathematics 2022,10, 3367 9 of 29 1 Ei =1 ˆ E 1 h4 i −νij Ei =−ˆ ν ˆ E 1 h2 ih2 j 1 Gij =1+ˆ ν ˆ E 1 h2 ih2 j , (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 modelAnisoH, we assumed ˆ ρ= 2.1 g/cm 3 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: H=ρ ˆ ρβ/4 A1/4 ˆ H1/2, (11) where ˆ H is the normalised fabric tensor. This tensor is normalised by imposing det(ˆ H) = 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 ρ/ˆ ρ 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: Y=2n2ˆ Gsym(HεH)(Hε)+ˆ λtr(H2ε)sym(Hε)o, (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: J=1−ω 3tr(Y)1+ωdev(Y), (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, Ψ∗ 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: gf(J,Ψ∗ t,w)≤0gr(J,Ψ∗ t,w)>0 resorption; gr(J,Ψ∗ t,w)≤0gf(J,Ψ∗ t,w)>0 formation; gf(J,Ψ∗ t,w)≤0gr(J,Ψ∗ t,w)≤0 lazy zone. (14) Mathematics 2022,10, 3367 16 of 29 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 Ψtis suitable. 3. Results The weighted Spearman correlation coefficients between bone density and the predictors proposed here are shown in Table 1for 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. 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 M S , and for the amplitude throughout the cycle, R A S . The predictors analysed are: the maximum principal stress ( σ1 ), the minimum principal stress ( σ3 ), the absolute maximum principal stress (AMPσ), the corresponding strains (ε1,ε3and AMPε), the strain energy density (SED or mechanical stimulus at the continuum level Ψ ), the mechanical stimulus at the tissue level ( Ψt ), the maximum tensile and compressive stresses ( σt and σc , respectively), the fluctuation of stresses ( σf ), the von Mises and Tresca stresses, the hydrostatic stress ( σo ) and the volumetric strain ( εo ). Values higher than 0.9 are highlighted in boldface; negative values are in red. Predictor GLT = 70 GLT = 80 GLT = 90 RM SRA SRM SRA SRM SRA S σ10.871 0.880 0.905 0.912 0.936 0.942 σ3−0.040 0.872 −0.047 0.897 −0.045 0.926 σt0.871 0.871 0.905 0.905 0.936 0.936 σc0.868 0.868 0.894 0.894 0.923 0.923 σf0.904 0.904 0.927 0.927 0.951 0.951 AMPσ0.896 0.896 0.921 0.920 0.945 0.944 σvonMises 0.898 0.898 0.922 0.922 0.948 0.948 σTresca 0.898 0.898 0.922 0.922 0.948 0.947 σo0.701 0.886 0.739 0.911 0.782 0.938 ε1−0.344 −0.340 −0.431 −0.429 −0.503 −0.502 ε30.410 −0.355 0.432 −0.445 0.457 −0.504 AMPε−0.328 −0.325 −0.427 −0.425 −0.499 −0.498 Ψ0.698 0.698 0.721 0.721 0.772 0.772 Ψtfor n=2 −0.048 −0.048 −0.092 −0.092 −0.097 −0.099 εo−0.202 −0.303 −0.255 −0.468 −0.327 −0.571 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, σM 3 , for obvious reasons, as the sign of σ3 (usually negative) is considered in the calculus of this peak. Therefore, σM 3 usually corresponds to the lowest absolute value throughout the cycle. The amplitude, σA 3 , is better correlated with the density as it usually measures the range of the compressive stress. Mathematics 2022,10, 3367 17 of 29 The strain magnitudes are poorly correlated with the density, in some cases with negative coefficients. SED (or Ψ ) is only moderately correlated with the density and not correlated at all if it is corrected to account for the porosity (Ψ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 2for 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 Afor 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 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 M P , and for the amplitude throughout the cycle, R A P . The predictors analysed in this case are: the absolute maximum principal stress (AMP σ ), the maximum tensile and compressive stresses ( σt and σc , respectively), the fluctuation of stresses ( σf ) and the mechanical stimulus at the continuum level ( Ψ ) and at the tissue level ( Ψ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 RM PRA PRM PRA PRM PRA P AMPσLinear 0.773 0.772 0.793 0.791 0.808 0.806 σtLinear 0.746 0.745 0.761 0.760 0.770 0.769 σcLinear 0.671 0.672 0.684 0.684 0.695 0.695 σfLinear 0.810 0.810 0.825 0.825 0.838 0.838 ΨLinear 0.012 0.012 0.012 0.012 0.012 0.012 Ψtfor n = 1 Linear −0.016 −0.016 −0.016 −0.016 −0.013 −0.013 Ψtfor n = 2 Linear −0.023 −0.023 −0.026 −0.026 −0.023 −0.023 AMPσQuadratic 0.784 0.782 0.804 0.803 0.818 0.817 σtQuadratic 0.762 0.761 0.778 0.778 0.787 0.786 σcQuadratic 0.678 0.679 0.690 0.691 0.700 0.700 σfQuadratic 0.823 0.823 0.839 0.839 0.850 0.850 ΨQuadratic 0.026 0.026 0.026 0.026 0.029 0.029 Ψtfor n = 1 Quadratic 0.017 0.017 0.017 0.017 0.018 0.018 Ψtfor n = 2 Quadratic 0.021 0.021 0.026 0.026 0.023 0.023 AMPσPower 0.780 0.778 0.799 0.798 0.813 0.811 σtPower 0.761 0.760 0.777 0.776 0.785 0.784 σcPower 0.674 0.675 0.683 0.684 0.690 0.691 σfPower 0.820 0.820 0.835 0.835 0.845 0.845 ΨPower 0.006 0.006 0.005 0.005 0.005 0.005 Ψtfor n = 1 Power −0.016 −0.016 −0.017 −0.017 −0.015 −0.015 Ψtfor n = 2 Power −0.019 −0.019 −0.022 −0.020 −0.015 −0.015 Tables 3and 4compare, 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) Mathematics 2022,10, 3367 18 of 29 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. 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 M S , and for the amplitude throughout the cycle, R A S . The predictors analysed are: the maximum principal stress ( σ1 ), the minimum principal stress ( σ3 ), the absolute maximum principal stress (AMPσ), the corresponding strains (ε1,ε3and AMPε), the strain energy density (SED or mechanical stimulus at the continuum level Ψ ), the mechanical stimulus at the tissue level ( Ψt ), the maximum tensile and compressive stresses ( σt and σc , respectively), the fluctuation of stresses ( σf ), the von Mises and Tresca stresses, the hydrostatic stress ( σo ) and the volumetric strain ( εo ). Values higher than 0.9 are highlighted in boldface; negative values are in red. Predictor AnisoH IsoH IsoJ RM SRA SRM SRA SRM SRA S σ10.905 0.912 0.906 0.913 0.905 0.913 σ3−0.047 0.897 −0.058 0.899 −0.059 0.895 σt0.905 0.905 0.906 0.906 0.905 0.905 σc0.894 0.894 0.895 0.895 0.892 0.892 σf0.927 0.927 0.928 0.928 0.925 0.925 AMPσ0.921 0.920 0.922 0.921 0.920 0.920 σvonMises 0.922 0.922 0.924 0.924 0.921 0.920 σTresca 0.922 0.922 0.923 0.923 0.920 0.920 σo0.739 0.911 0.742 0.913 0.739 0.911 ε1−0.431 −0.429 −0.431 −0.429 −0.443 −0.441 ε30.432 −0.445 0.431 −0.443 0.429 −0.458 AMPε−0.427 −0.425 −0.425 −0.423 −0.440 −0.438 Ψ0.721 0.721 0.722 0.721 0.708 0.708 Ψtfor n = 2 −0.092 −0.092 −0.092 −0.092 −0.114 −0.114 εo−0.255 −0.468 −0.254 −0.464 −0.247 −0.489 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). 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 ( σf , AMP σ ) and only the coefficients for the maximum variable throughout the cycle, R M S . The other stress variables follow the same trend, as well as R A S . 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. Mathematics 2022,10, 3367 19 of 29 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 M P , and for the amplitude throughout the cycle, R A P . The predictors analysed in this case are: the absolute maximum principal stress (AMPσ), the maximum tensile and compressive stresses (σtand σc , respectively), the fluctuation of stresses ( σf ), the mechanical stimulus at the continuum level ( Ψ ) and at the tissue level ( Ψt ) for two values of n . Values higher than 0.8 are highlighted in boldface; negative values are in red. Predictor Type of AnisoH IsoH IsoJ Correlation RM PRA PRM PRA PRM PRA P AMPσLinear 0.793 0.791 0.802 0.800 0.814 0.813 σtLinear 0.761 0.760 0.767 0.766 0.771 0.770 σcLinear 0.684 0.684 0.692 0.692 0.707 0.708 σfLinear 0.825 0.825 0.832 0.832 0.840 0.840 ΨLinear 0.012 0.012 0.012 0.012 0.015 0.015 Ψtfor n = 1 Linear −0.016 −0.016 −0.016 −0.016 −0.017 −0.017 Ψtfor n = 2 Linear −0.026 −0.026 −0.026 −0.026 −0.029 −0.029 AMPσQuadratic 0.804 0.803 0.814 0.812 0.834 0.833 σtQuadratic 0.778 0.778 0.785 0.784 0.797 0.796 σcQuadratic 0.690 0.691 0.698 0.699 0.719 0.720 σfQuadratic 0.839 0.839 0.846 0.846 0.861 0.861 ΨQuadratic 0.026 0.026 0.026 0.026 0.030 0.030 Ψtfor n = 1 Quadratic 0.017 0.017 0.017 0.017 0.020 0.020 Ψtfor n = 2 Quadratic 0.026 0.026 0.026 0.026 0.029 0.029 AMPσPower 0.799 0.798 0.809 0.808 0.833 0.832 σtPower 0.777 0.776 0.783 0.783 0.797 0.796 σcPower 0.683 0.684 0.692 0.693 0.717 0.718 σfPower 0.835 0.835 0.842 0.842 0.861 0.861 ΨPower 0.005 0.005 0.006 0.006 0.008 0.008 Ψtfor n = 1 Power −0.017 −0.017 −0.017 −0.017 −0.018 −0.018 Ψtfor n = 2 Power −0.022 −0.022 −0.022 −0.022 −0.023 −0.023 Table 5. Weighted Pearson correlation coefficients (power fit) obtained for the maximum variable throughout the cycle, R M S , for GLT = 80, constitutive model IsoH and using C3D4 elements. The predictors analysed are: the fluctuation of stresses ( σf ) and the absolute maximum principal stress (AMP σ ). Values higher than 0.8 are highlighted in boldface. Proximal Diaphysis Distal Global σf0.718 0.879 0.660 0.844 AMPσ0.694 0.862 0.640 0.813 The influence of the mesh (C3D4 vs. C3D10) was analysed by comparing the correlation coefficients in Table 6; 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 Ψ 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. Mathematics 2022,10, 3367 20 of 29 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 M P , and for the amplitude throughout the cycle, R A P . 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 C3D4 C3D10 Correlation RM PRA PRM PRA P σ1Spearman 0.906 0.913 0.904 0.912 σ3Spearman −0.058 0.899 −0.236 0.896 σtSpearman 0.906 0.906 0.904 0.904 σcSpearman 0.895 0.895 0.893 0.894 σfSpearman 0.928 0.928 0.930 0.930 AMPσSpearman 0.922 0.921 0.923 0.922 σvonMises Spearman 0.924 0.924 0.925 0.925 σTresca Spearman 0.923 0.923 0.925 0.924 σoSpearman 0.742 0.913 0.734 0.911 ε1Spearman −0.431 −0.429 −0.499 −0.496 ε3Spearman 0.431 −0.443 0.484 −0.495 AMPεSpearman −0.425 −0.423 −0.482 −0.480 ΨSpearman 0.722 0.721 0.709 0.709 Ψtfor n = 2 Spearman −0.092 −0.092 −0.182 −0.182 εoSpearman −0.254 −0.464 −0.287 −0.499 AMPσPearson (Power) 0.809 0.808 0.874 0.874 σtPearson (Power) 0.783 0.783 0.818 0.818 σcPearson (Power) 0.692 0.693 0.741 0.742 σfPearson (Power) 0.842 0.842 0.902 0.902 ΨPearson (Power) 0.006 0.006 0.011 0.011 Ψtfor n = 1 Pearson (Power) −0.017 −0.017 −0.019 −0.019 Ψtfor n = 2 Pearson (Power) −0.022 −0.022 −0.025 −0.025 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 σ or AMP ε are within a given range. The elements of cortical ( ρ> 1.2 g/cm 3 ) and trabecular bone ( ρ≤ 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 ε is in the range [200,600] µε occupy 30% of the total volume of trabecular bone (ρ≤1.2 g/cm3). 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 ε in the range of 200–600 µε . This is not so evident in trabecular bone, though 51% is still within the range of 200–1000 µε and 65% is in the range of 200–1400 µε . This is still a narrow range, as overload strains are up to 4000 µε[36,37]. Mathematics 2022,10, 3367 21 of 29 AMPεtrabecular 0 5 10 15 20 25 30 35 % µ 200–600 600–1000 1000–1400 1400–1800 1800–2200 >2200 0–200 AMPεcortical 200–600 0 10 20 30 40 50 60 70 80 90 600–1000 1000–1400 1400–1800 1800–2200 >2200 % µε 0–200 AMPσtrabecular 0 5 10 15 20 25 30 35 40 45 50 0–0.001 0.001–0.01 0.01–0.1 0.1–1 1–2 >2 % MPa AMPσcortical 0 5 10 15 20 25 30 35 40 45 50 0–2 2–4 4–6 6–8 8–10 10–12 >12 % MPa IsoH IsoJ AnisoH Figure 6. Distribution of the absolute maximum principal stress (AMP σ ) and strain (AMP ε ) for the elements of trabecular ( ρ≤ 1.2 g/cm 3 ) and cortical bone ( ρ< 1.2 g/cm 3 ). The Y axis shows the percentage of the total volume occupied by the elements, whose AMP σ or AMP ε are in a given range. 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 – 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 Tables 1and 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 σ ) and the fluctuation of stresses ( σ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 ρ . 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é ( Ψt ) is even worse, as it yields negative correlation coefficients. Mathematics 2022,10, 3367 22 of 29 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 σf and AMP σ showed a value of R P slightly over 0.8 (up to 0.861 for σf in the IsoJ model). This means that R 2>∼ 0.64, i.e., around 64% of the variation of the density, can be explained by σf(up to 74% in the IsoJ model). Only some of the variables appearing in Table 1were chosen for assessing its Pearson coefficient, those having a high Spearman coefficient plus SED and its related variables ( Ψ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 σ ) and one representing the strain-state (AMP ε ), was analysed by means of histograms (see Figure 6) that distinguish between cortical ( ρ> 1.2 g/cm 3 ) and trabecular bone ( ρ≤ 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 ε∈ [200–1400] µε 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ε∈[200–600] µε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 µε [ 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 µε [ 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−ρ , is being forced by the Mathematics 2022,10, 3367 23 of 29 relationship between stiffness and density, E=F(ρ) , in particular for those predictors derived from the stress tensor (for the sake of generality, we should write C=F(ρ) , with 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(ρ) ) . 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 σ=H(E)to be monotonically increasing as well. Thus, if we write: σ=HF(ρ), (38) we can expect a monotonically increasing function relating σ and ρ , in other words, a positive correlation S−ρ for the predictors derived from the stress tensor. Therefore, it could seem that the function E=F(ρ) is forcing S−ρ and it is so to some extent, especially for the Spearman correlation coefficients. However, if E=F(ρ) were the only determinants of the correlation S−ρ , 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 σ histograms, but both are almost identical for the three models (see Table 3and 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−ρ 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: σ=HF(ρ),F, (39) where 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(σ), thus leading to: S=GHF(ρ),F. (40) In view of Equation (40) the predictors, S , need not be predetermined only by the function F(ρ) , 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 Tables 1and 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−ρ . 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 Si e (recall Section 2.7) and these variables were correlated with ρe to give R i P , a weighted Pearson correlation coefficient for each instant i . The worst coefficient between AMP σ and density throughout the cycle was R i P = 0.48 (for the quadratic fit and in the case GLT = 80, constitutive model IsoH and C3D4 elements), which is far from R M P = 0.804. The former is a very poor correlation. (R i P)2 = 0.23, i.e., only about 25% of the variance of density can be explained by AMP σ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. Mathematics 2022,10, 3367 24 of 29 In summary, the function E=F(ρ) is contributing to obtain the correlations S−ρ , 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. Mathematics 2022,10, 3367 25 of 29 Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/math10183367/s1. Author Contributions: Conceptualization, J.M.-R., J.O., J.L.C.-G. and M.Á.P.; methodology, J.M.-R., J.O. and J.L.C.-G.; software, F.G.-M., J.L.C.-G., J.M.-R. and J.O.; formal analysis, F.G.-M., J.M.-R. and J.O.; writing—original draft preparation, J.M.-R., J.O. and J.L.C.-G; writing—review and editing, J.L.C.-G. and M.Á.P.; visualization, J.L.C.-G., J.O. and J.M.-R.; funding acquisition, J.M.-R. All authors have read and agreed to the published version of the manuscript. Funding: Funding was provided by the Fondo Europeo de Desarrollo Regional (FEDER), the Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía, dentro del Programa Operativo FEDER 2014–2020 and the Universidad de Sevilla to the project P18-RT-3611 entitled “Efecto combinado del ejercicio físico y el denosumab en el tratamiento de la osteoporosis. Diseño de un tratamiento farmacológico específico de paciente” for which this paper was prepared. Institutional Review Board Statement: The study was conducted in accordance with the Declaration of Helsinki, and approved by the Andalusian Biomedical Research Ethics Platform (approval number 20151012181252). Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Data Availability Statement: The finite element model, including boundary conditions and loads, used in this study is available in the Supplementary Material. Other data presented in this study are not publicly available because they contain patient information. Conflicts of Interest: The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Abbreviations The following abbreviations are used in this manuscript: BRM Bone remodelling model(s) SED Strain energy density FE(M) Finite Element (Method) CT Computer tomography DPTBR Density prediction through bone remodelling GLT Grey level threshold AMPεAbsolute maximum principal strain AMPσAbsolute maximum principal stress MF Muscle force(s) HJF Hip joint force KJF Knee joint force Appendix A. Weighted Pearson and Spearman Correlation Coefficients In the case of a linear correlation, the Pearson coefficient between the bone apparent density ρ and the predictor var is weighted with the volume of the finite elements, which are the sample points. The weighted Pearson coefficient is defined as the weighted covariance divided by the weighted variances of both variables and determines the proportion of variance explained by the weighted linear fit of the point cloud (ρe,vare): RP= n ∑ e=1 [Ve(ρe−ρ)(vare−var)] sn ∑ e=1 [Ve(ρe−ρ)2]· n ∑ e=1 [Ve(vare−var)2] (A1)