**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 1 and 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 (*<sup>σ</sup>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.

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<sup>2</sup> > ∼ 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 1 were 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*<sup>S</sup>*. 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/cm3) and trabecular bone (*ρ* ≤ 1.2 g/cm3). 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 sugges<sup>t</sup> 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 ye<sup>t</sup> 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 rela-

tionship between stiffness and density, *E* = *<sup>F</sup>*(*ρ*), in particular for those predictors derived from the stress tensor (for the sake of generality, we should write **C** = *<sup>F</sup>*(*ρ*), 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* = *<sup>F</sup>*(*ρ*)) . 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:

$$
\sigma = H(F(\rho)),
\tag{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* = *<sup>F</sup>*(*ρ*) is forcing *S* − *ρ* and it is so to some extent, especially for the Spearman correlation coefficients. However, if *E* = *<sup>F</sup>*(*ρ*) 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 ye<sup>t</sup> 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 3 and Figure 6). Hence, we can conclude that the assumed constitutive model has no significant effect on the correlations.

However, there is another argumen<sup>t</sup> 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:

$$
\sigma = H(F(\rho), \mathbf{F}),
\tag{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* = *<sup>G</sup>*(*σ*), thus leading to:

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

In view of Equation (40) the predictors, *S*, need not be predetermined only by the function *<sup>F</sup>*(*ρ*), 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 1 and 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.

In summary, the function *E* = *<sup>F</sup>*(*ρ*) 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.
