Bone remodelling models (BRM) are often used to estimate the density distribution in bones from the loads they are subjected to. BRM deﬁne 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; ﬂuctuation 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, ˙r: ˙r=f(S) Evolutionary model (1) where Sis 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 deﬁned 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 signiﬁcant. 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 ﬁnal 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 Sand 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. r=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 evolution- ary 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 Sis the variable that drives bone remodelling and must comprise bi- ological 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 inﬂuence 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 deﬁned 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 brieﬂy 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 brieﬂy describe how the gait cycle was simulated with a FE model. The correlation coefﬁcients used to evaluate the relationship between bone density and the variables used as candidates for mechanical stimuli (or predictors) are deﬁned in Section 2.8. In Section 2.9, we deﬁne 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, deﬁned 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 Depart- ment of Mechanical Engineering and Manufacturing of the Universidad de Sevilla. A Vi- con®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 modiﬁed Cleveland proto- col [16]. The subject under study was an adult male of 27 years old, with no reported patholo- gies, 1.85 m tall and weighing 75 kg, who walked at a freely chosen forward speed. The par- ticipant 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 01234 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 ﬁt 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 labora- tory. 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 dy- namics 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. Addition- ally, 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 ﬁnal 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]: re=A+Bge, (3) where geandreare the greyscale value and the density estimated for element e, respectively. Two pairs of greyscale–density values are needed to deﬁne the linear relationship (constants Aand 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 ﬁrst 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) thatre=0forge=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: re=2.1geGLT 255GLT. (4) The greyscale distribution was mapped from the CT scan to the FE mesh using the software Bonemat , and Equation (3)was used to convert the greyscale value into bone apparent density. Those elements with a greyscale value below GLT were assigned a very small density (0.001 g/cm3), thus, yielding a negligible—although, not null—stiffness, something necessary to avoid convergence problems. Figure 2 shows the distributions of density obtained for the three GLT. It can be seen that those three GLT led to distributions of the estimated density which are similar, in general, with the exception of the proximal region and the thickness of the cortical layer in the diaphysis. As GLT rises, the volume identiﬁed 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 deﬁned 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 n= 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) =2014r2.5ifr1.2 g/cm3 1763r3.2ifr>1.2 g/cm3 (5a) n=0.2 if r1.2 g/cm3 0.32 if r>1.2 g/cm3. (5b) Hernandez et al. [ 29] proposed another correlation for E, based on the ash fraction, a(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=r ˆr, with rand ˆrbeing the apparent density and the tissue density, respectively, then the correlation proposed by Hernandez et al. reads: E(MPa) =84370r ˆr2.58 a2.74. (6) If typical values of ash fraction a=0.68 and tissue density ˆr=2.31 g/cm3[29] are used, Equation (6) becomes: EHernandez (MPa) =3388r2.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 n=0.3was used in conjunction with Equation (7) in model IsoH . 0 0.5 1 1.5 2 2.50510152025 0 0.5 1 1.5 2 2.500.511.52 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 proper- ties 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 adap- tation 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 deﬁning a stepwise Mathematics 2022 ,10, 3367 8 of 29 linear relationship between the mechanical stimulus, Yt, and the bone formation/resorption rate, ˙r(see Figure 4). The reference value of the stimulus, Y tdeﬁnes 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 ˙rdetermines the temporal evolution of apparent density, ˙r. In turn, apparent density determines the variation of the Young’s modulus through correlations such as (5a) or (7). Ψt−Ψ∗ t˙r 1c2 1 c1w 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, Yt, and a reference value of that stimulus, Y 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= (1D)C0, (8) where Dis 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 Dis replaced with a damage tensor Dand 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=1H2, (9) with 1being the second order identity tensor. Equation (9)leads to the following relation- ship 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 ˆE1 h4 i nij Ei=ˆn ˆE1 h2 ih2 j 1 Gij=1+ˆn ˆE1 h2 ih2 j,(10) where ˆEand ˆnare, 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 requal to the density of the bone matrix, ˆr, which is assumed constant. In this particular model AnisoH , we assumed ˆr=2.1g/cm3and chose the Hernandez correlation (7), together with n=0.3. The inﬂuence of porosity and directionality are factorised in this model by redeﬁning the fabric tensor as follows: H=r ˆrb/4 A1/4ˆH1/2, (11) where ˆHis 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, bis the exponent of the apparent density in the correlations (5a) or(7)and Ais 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. Acan here be considered a constant. As stated before, the quotient r/ˆrin Equation (11) is equal to the bone volume fraction, BV/TV=1p, with pas the porosity. The mechanical stimulus is deﬁned in this model through the tensor, Y: Y=2n 2ˆGsym(H#H)(H#)+ˆltr(H2#)sym(H#)o , (12) where ˆGand ˆlare 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 deﬁned another tensor, J, to quantify the relative inﬂuence of the spherical and deviatoric parts of the stimulus: J=1w 3tr(Y)1+wdev(Y), (13) where 1is the identity tensor, dev()represents the deviatoric part of a tensor and the anisotropy factor, w, weights the importance of the anisotropy of the stimulus in the model. This factor ranges from w= 0, which means that the model only depends on the isotropic component of the stimulus, to w= 1, which produces the maximum level of anisotropy. The same value used by Doblaré and García [ 4] was used here ( w= 0.1). Two functions grand gfare proposed to establish the remodelling criteria. These functions depend on the stimulus Jand 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, Y 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,Y t, w)0gr(J,Y t, w)>0 resorption; gr(J,Y t, w)0gf(J,Y t, w)>0 formation; gf(J,Y t, w)0gr(J,Y t, w)0 lazy zone.(14) Mathematics 2022 ,10, 3367 10 of 29 Based on the fulﬁlment of the corresponding criterion, the variation of the fabric tensor H, that accounts for the variation of anisotropy (through tensor ˆH) and porosity (see Equation (11)), is provided by: ˙H=k1ˆr rJ3ˆwresorption; ˙H=k2ˆr rJ ˆw formation; ˙H=0 lazy zone;(15) where the tensor ˆwis introduced to simplify the expression, as follows: ˆw=12w 31 1+wI (16) with Ibeing the fourth-order identity tensor. The factors k1and k2in Equation (15) de- pended on several parameters in the original model. One of those parameters is ˙r, so that the amount of formed or resorbed tissue modiﬁes the fabric tensor through porosity. 2.5. Modiﬁcation of the BRM to Maintain Density Constant In the original model, Doblaré and García used ˙Hto 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 k1and k2can be assumed as constants. In such case, by deriving Equation (11): ˙H=1 2r ˆrb/4 A1/4ˆH1/2˙ˆH (17) and given that ˆHmust remain normalised (det (ˆH) =1), we ﬁnally adopted: ˙ˆH=cˆH1/2˙H, (18) where cis a constant. This expression can be used in an Euler forward integration algorithm to yield: ˆH(tj+1) = ˆH(tj) +cˆH1/2(tj)˙H(tj)Dt, (19) where tj+1and tjare two consecutive integration steps, the time step Dt=1is chosen and cis the constant necessary to enforce the condition det ˆH(tj+1) =1. The simulation started from an initially isotropic material ( ˆH(t1)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.: ån e=1ˆHe(tj+1):ˆHe(tj+1)ån e=1ˆHe(tj):ˆHe(tj) ån e=1ˆHe(tj):ˆHe(tj)<0.001. (20) 2.6. Modiﬁcation 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 modiﬁed 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 modiﬁcation. 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 Mathematics 2022 ,10, 3367 11 of 29 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. Ψ∗ t tFormation A1 gf<0gf>0 ResorptionA3gr>0gr<0Lazy zoneA2Peaks:gM fandgm r Valleys:gm fandgM r Figure 5. Remodelling criteria with different loads. The horizontal lines limit the zones of formation (top line) and resorption (bottom line). If the peak is over the top line, then gM f>0and gmr<0 and formation occurs. If the peak is below the bottom line, then gM f<0and gmr>0and resorption occurs. If the peak lies between both lines, gM f<0and gmr<0and neither formation nor resorption occurs (lazy zone). The peaks of the stimulus coincide with the maximum of the formation criterion function, gf, which is proportional to the stimulus. Those peaks are termed here, gM f. Since the resorption criterion function, gr, is inversely proportional to the stimulus [4], the local minimum of this function, gm r, also coincides with the peaks of the stimulus (it must be noted that gM fand gm rdo not necessarily coincide. They are simply reached at the same instant). Analogously, the valleys of the stimulus coincide with the minimum of gfand the maximum of gr, termed here gm fand gM r, 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 gfand gr) at every point of the FE mesh throughout the cycle, in order to capture the peaks gM fand gm rfor 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. Mathematics 2022 ,10, 3367 12 of 29 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: gM f>0 formation; gm r0gM f0 lazy zone; gM f0gm r>0 resorption.(21) It is important to highlight that gM fand gm rare 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 modiﬁcation 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 simpliﬁcation 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 WR3be an open bounded domain and G= ¶Wbe its boundary, assumed to be Lipschitz continuous and divided into two disjoint parts GDandGNwhere Dirichlet and Neumann boundary conditions are applied, respectively. We denote, by x=xiei, a generic point of Wand n(x) =nieias the outward unit normal vector to Gat a point x. The Einstein summation notation was adopted and a Cartesian basis ei(i=1, 2, 3) can be used without loss of generality. Letu=uiei,sss=sijei ejand###(u) =#ijei ejdenote the displacement ﬁeld, the stress tensor and the linearised strain tensor, respectively. ei ej(i,j=1, 2, 3 ) represents the Cartesian tensorial basis. Let b=bieidenote the known vector of body forces, t=tiei the known vector of surface traction forces at GNandddd=dieithe known displacements atGD. The Momentum Conservation Principle states: sij,j+bi=r¨ui, (22) where ,jdenotes the partial derivative ¶/¶xjand¨()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: #ij(u) =1 2 ui,j+uj,i i,j=1, 2, 3 (23) and ﬁnally, strains and stresses are related by the constitutive equation, which for linear elastic materials reads: sij=Cijkl#kl i,j,k,l=1, 2, 3, (24) where Cijklare the components of the fourth-rank stiffness tensor C. In the case of an isotropic material, this tensor is completely deﬁned by the Young’s modulus, E, and the Poisson’s ratio, n, 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 GDandGN: u=dddatGD (25a) sijnj=tiatGN. (25b) Mathematics 2022 ,10, 3367 13 of 29 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 deﬁnition 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 deﬁned 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 deﬁned as passing through the corresponding joint centre (deﬁned 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 ( r¨uiin 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 Material. This pseudo-static analysis provided the temporal evolution of stresses and strains at each element eof 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, identiﬁed with the element. In the case of C3D10 elements, the variables were evaluated at the centroid of the element) and at every instant iof 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 iand for each element, e, thus, yielding Si e. Then, the following maximum, minimum and amplitude were deﬁned to represent its evolution throughout the gait cycle: SM e=Max i(Si e) (26a) Sm e=Min i(Si e) (26b) SA e=SM eSm e. (26c) Note that the deﬁnition of SM eis 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 modiﬁcation of the BRM presented in Section 2.6 is related to this idea. The variable, SM e, generalises this concept of the peak of stress to the peak of stimulus. As an alternative to the maximum of the stimulus throughout the cycle, the amplitude is considered in SA e. 2.8. Correlation Coefﬁcients The Spearman and Pearson correlation coefﬁcients were used to assess the statistical dependence between the bone density and the variables deﬁned later in Section 2.9. Based on the deﬁnitions made in Equation (26), the following coefﬁcients were calculated, taking each element as a point of the sample: • RM j, between the maximum throughout the cycle SM eand the apparent density re. Mathematics 2022 ,10, 3367 14 of 29 • RA j, between the amplitude throughout the cycle SA eand the apparent density re. where jstands for P(Pearson) or S(Spearman). Weighted coefﬁcients 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 coefﬁcients are given in Appendix A. The Spearman correlation coefﬁcient 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 coefﬁcient is a parametric measure of correlation between two variables that assesses if they are related by a speciﬁc function. In particular, the linear, quadratic and power correlation coefﬁcients were calculated for those variables that yielded a high Spearman coefﬁcient, in order to conﬁrm 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 coefﬁcient, it was also correlated with bone density through the Pearson coefﬁcient 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—Deﬁnition 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, s1s2s3, and analogously, #1#2#3. The maximum and minimum principal stresses ( s1,s3) and strains ( #1,#3) are proposed as predictor variables. The maximum tensile stress is deﬁned as: st=s1ifs10 0 if s1<0(27) and the maximum compressive stress as: sc=s3ifs3<0 0 if s30. (28) The ﬂuctuations of stresses throughout the cycle can be measured by the variable: sf=sM 1sm 3, (29) where the superscripts Mand mfollow the deﬁnition given in Equation (26). The absolute maximum principal stress (AMP s) and strain (AMP #) are deﬁned analogously as: AMP s=Maxjs1j,js3j (30a) AMP #=Maxj#1j,j#3j . (30b) The von Mises ( svonMises ) and Tresca ( sTresca ) 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 ( so) and volumetric strain ( #o) are also proposed as predictor variables: svonMises =r (s1s2)2+ (s1s3)2+ (s2s3)2 2(31a) sTresca =s1s3 2(31b) so=s1+s2+s3 3(31c) Mathematics 2022 ,10, 3367 15 of 29 #o=#1+#2+#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 ( sss) and strain ( ###) tensors or their components ( sij,#ij): SED=1 2sss:###=1 2sij#ij. (32) Beaupré et al. [ 3] deﬁned the mechanical stimulus that drives bone remodelling, Y, 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 Yand SED [3]: Y=kSED. (33) In other words, we can identify Yand SED for correlation purposes. More importantly, Beaupré et al. proposed to take into account the porosity of the tissue to redeﬁne the mechanical stimulus. Thus, Yrepresents the mechanical stimulus at the continuum (or macroscopic) level. SED can be calculated through FEM and if the constant kin Equation (33) is known, the mechanical stimulus at the continuum level, Y, 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, Ycan 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]: Yt=Y (1p)n, (34) where p2[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 coefﬁcients 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 sss=seeei eeei, with eeeibeing the loading direction and sthe applied stress. In such case, and assuming a linearly elastic and isotropic behaviour for bone, the SED can be calculated through Equation (32) as: SED=1 2s #=1 2E#2, (35) where ###is the strain tensor and #is the strain in the loading direction. The asterisk has been added to highlight that this expression corresponds to a particular case. Additionally, a typical power correlation between the Young’s modulus and the apparent bone density can be assumed, for example Equations (5a) and (7), which would read: E=Brb, (36) where Bandbare constants. In this case, Ycan be rewritten using Equations (33) and (35) as : Y=k1 2B#2rb=Krb, (37) Mathematics 2022 ,10, 3367 16 of 29 where the constants preceding the factor rbwere 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, Y could be related to the bone apparent density through a power law. Recalling Equation (34) and given that (1p)is equal to the bone volume fraction, which is proportional to the bone density, the mechanical stimulus at the tissue level, Yt, 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 Y(or equivalently SED) or Ytis suitable. 3. Results The weighted Spearman correlation coefﬁcients between bone density and the predic- tors proposed here are shown in Table 1 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 coefﬁcients are given in the three cases for comparison. The fact that the Spearman correlation coefﬁcient 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 coefﬁcients obtained with the constitutive model, AnisoH , using C3D4 elements and for different values of GLT: for the maximum variable throughout the cycle, RM S, and for the amplitude throughout the cycle, RA S. The predictors analysed are: the maximum principal stress ( s1), the minimum principal stress ( s3), the absolute maximum principal stress (AMP s), the corresponding strains ( #1,#3and AMP #), the strain energy density (SED or mechanical stimulus at the continuum level Y), the mechanical stimulus at the tissue level ( Yt), the maximum tensile and compressive stresses ( standsc, respectively), the ﬂuctuation of stresses ( sf), the von Mises and Tresca stresses, the hydrostatic stress ( so) and the volumetric strain ( #o). Values higher than 0.9 are highlighted in boldface; negative values are in red. PredictorGLT = 70 GLT = 80 GLT = 90 RM S RA S RM S RA S RM S RA S s1 0.871 0.880 0.905 0.912 0.936 0.942 s30.040 0.872 0.047 0.897 0.045 0.926 st 0.871 0.871 0.905 0.905 0.936 0.936 sc 0.868 0.868 0.894 0.894 0.923 0.923 sf 0.904 0.904 0.927 0.927 0.951 0.951 AMP s 0.896 0.896 0.921 0.920 0.945 0.944 svonMises 0.898 0.898 0.922 0.922 0.948 0.948 sTresca 0.898 0.898 0.922 0.922 0.948 0.947 so 0.701 0.886 0.739 0.911 0.782 0.938 #10.3440.3400.4310.4290.5030.502 #3 0.410 0.355 0.432 0.445 0.457 0.504 AMP #0.3280.3250.4270.4250.4990.498 Y 0.698 0.698 0.721 0.721 0.772 0.772 Ytfor n=20.0480.0480.0920.0920.0970.099 #o0.2020.3030.2550.4680.3270.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, sM 3, for obvious reasons, as the sign ofs3(usually negative) is considered in the calculus of this peak. Therefore, sM 3usually corresponds to the lowest absolute value throughout the cycle. The amplitude, sA 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 coefﬁcients. SED (or Y) is only moderately correlated with the density and not correlated at all if it is corrected to account for the porosity ( Yt). Regarding the inﬂuence 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 coefﬁcients are shown in Table 2 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 coefﬁcient together with the SED at the continuum and at the tissue level for n = 1 and n = 2 (see Equation (34)). The Pearson coefﬁcients 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 ﬁt 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 coefﬁcients were obtained, as in the case of the Spearman coefﬁcients. Table 2. Weighted Pearson correlation coefﬁcients obtained with the constitutive model, AnisoH , using C3D4 elements and for different values of GLT: for the maximum variable throughout the cycle, RM P, and for the amplitude throughout the cycle, RA P. The predictors analysed in this case are: the absolute maximum principal stress (AMP s), the maximum tensile and compressive stresses (standsc, respectively), the ﬂuctuation of stresses ( sf) and the mechanical stimulus at the continuum level ( Y) and at the tissue level ( Yt) for two values of n. Values higher than 0.8 are highlighted in boldface; negative values are in red. Predictor TypeGLT = 70 GLT = 80 GLT = 90 RM P RA P RM P RA P RM P RA P AMP s Linear 0.773 0.772 0.793 0.791 0.808 0.806 st Linear 0.746 0.745 0.761 0.760 0.770 0.769 sc Linear 0.671 0.672 0.684 0.684 0.695 0.695 sf Linear 0.810 0.810 0.825 0.825 0.838 0.838 Y Linear 0.012 0.012 0.012 0.012 0.012 0.012 Ytfor n = 1 Linear 0.0160.016 0.0160.016 0.0130.013 Ytfor n = 2 Linear 0.0230.023 0.0260.026 0.0230.023 AMP s Quadratic 0.784 0.782 0.804 0.803 0.818 0.817 st Quadratic 0.762 0.761 0.778 0.778 0.787 0.786 sc Quadratic 0.678 0.679 0.690 0.691 0.700 0.700 sf Quadratic 0.823 0.823 0.839 0.839 0.850 0.850 Y Quadratic 0.026 0.026 0.026 0.026 0.029 0.029 Ytfor n = 1 Quadratic 0.017 0.017 0.017 0.017 0.018 0.018 Ytfor n = 2 Quadratic 0.021 0.021 0.026 0.026 0.023 0.023 AMP s Power 0.780 0.778 0.799 0.798 0.813 0.811 st Power 0.761 0.760 0.777 0.776 0.785 0.784 sc Power 0.674 0.675 0.683 0.684 0.690 0.691 sf Power 0.820 0.820 0.835 0.835 0.845 0.845 Y Power 0.006 0.006 0.005 0.005 0.005 0.005 Ytfor n = 1 Power 0.0160.016 0.0170.017 0.0150.015 Ytfor n = 2 Power 0.0190.019 0.0220.020 0.0150.015 Tables 3 and 4 compare, respectively, the Spearman and Pearson coefﬁcients obtained using the three constitutive models analysed in this work: AnisoH ,IsoH and IsoJ. As in- dicated 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 ﬁt 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 coefﬁcients obtained for GLT = 80, constitutive models AnisoH ,IsoH and IsoJand using C3D4 elements: for the maximum variable throughout the cycle, RM S, and for the amplitude throughout the cycle, RA S. The predictors analysed are: the maximum principal stress ( s1), the minimum principal stress ( s3), the absolute maximum principal stress (AMP s), the corresponding strains ( #1,#3and AMP #), the strain energy density (SED or mechanical stimulus at the continuum level Y), the mechanical stimulus at the tissue level ( Yt), the maximum tensile and compressive stresses ( standsc, respectively), the ﬂuctuation of stresses ( sf), the von Mises and Tresca stresses, the hydrostatic stress ( so) and the volumetric strain ( #o). Values higher than 0.9 are highlighted in boldface; negative values are in red. PredictorAnisoH IsoH IsoJ RM S RA S RM S RA S RM S RA S s1 0.905 0.912 0.906 0.913 0.905 0.913 s30.047 0.897 0.058 0.899 0.059 0.895 st 0.905 0.905 0.906 0.906 0.905 0.905 sc 0.894 0.894 0.895 0.895 0.892 0.892 sf 0.927 0.927 0.928 0.928 0.925 0.925 AMP s 0.921 0.920 0.922 0.921 0.920 0.920 svonMises 0.922 0.922 0.924 0.924 0.921 0.920 sTresca 0.922 0.922 0.923 0.923 0.920 0.920 so 0.739 0.911 0.742 0.913 0.739 0.911 #10.4310.429 0.4310.429 0.4430.441 #3 0.4320.445 0.431 0.443 0.429 0.458 AMP #0.4270.425 0.4250.423 0.4400.438 Y 0.721 0.721 0.722 0.721 0.708 0.708 Ytfor n = 20.0920.092 0.0920.092 0.1140.114 #o0.2550.468 0.2540.464 0.2470.489 We have also analysed the spatial distribution of the correlations, in particular of the Pearson coefﬁcients (power ﬁt), 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 inﬂuence 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 (sf, AMP s) and only the coefﬁcients for the maximum variable throughout the cycle, RM S. The other stress variables follow the same trend, as well as RA 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 coefﬁcients are high in the diaphysis, signiﬁcantly worse in the proximal and especially worse in the distal third, inﬂuenced by the simpliﬁed 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 simpliﬁcation 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 coefﬁcients obtained for GLT = 80, constitutive models AnisoH, IsoH and IsoJand using C3D4 elements: for the maximum variable throughout the cycle, RM P, and for the amplitude throughout the cycle, RA P. The predictors analysed in this case are: the absolute maximum principal stress (AMP s), the maximum tensile and compressive stresses ( stand sc, respectively), the ﬂuctuation of stresses ( sf), the mechanical stimulus at the continuum level ( Y) and at the tissue level ( Yt) for two values of n. Values higher than 0.8 are highlighted in boldface; negative values are in red. PredictorType of AnisoH IsoH IsoJ Correlation RM P RA P RM P RA P RM P RA P AMP s Linear 0.793 0.791 0.802 0.800 0.814 0.813 st Linear 0.761 0.760 0.767 0.766 0.771 0.770 sc Linear 0.684 0.684 0.692 0.692 0.707 0.708 sf Linear 0.825 0.825 0.832 0.832 0.840 0.840 Y Linear 0.012 0.012 0.012 0.012 0.015 0.015 Ytfor n = 1 Linear 0.0160.016 0.0160.016 0.0170.017 Ytfor n = 2 Linear 0.0260.026 0.0260.026 0.0290.029 AMP s Quadratic 0.804 0.803 0.814 0.812 0.834 0.833 st Quadratic 0.778 0.778 0.785 0.784 0.797 0.796 sc Quadratic 0.690 0.691 0.698 0.699 0.719 0.720 sf Quadratic 0.839 0.839 0.846 0.846 0.861 0.861 Y Quadratic 0.026 0.026 0.026 0.026 0.030 0.030 Ytfor n = 1 Quadratic 0.017 0.017 0.017 0.017 0.020 0.020 Ytfor n = 2 Quadratic 0.026 0.026 0.026 0.026 0.029 0.029 AMP s Power 0.799 0.798 0.809 0.808 0.833 0.832 st Power 0.777 0.776 0.783 0.783 0.797 0.796 sc Power 0.683 0.684 0.692 0.693 0.717 0.718 sf Power 0.835 0.835 0.842 0.842 0.861 0.861 Y Power 0.005 0.005 0.006 0.006 0.008 0.008 Ytfor n = 1 Power 0.0170.017 0.0170.017 0.0180.018 Ytfor n = 2 Power 0.0220.022 0.0220.022 0.0230.023 Table 5. Weighted Pearson correlation coefﬁcients (power ﬁt) obtained for the maximum variable throughout the cycle, RM S, for GLT = 80, constitutive model IsoH and using C3D4 elements. The predic- tors analysed are: the ﬂuctuation of stresses ( sf) and the absolute maximum principal stress (AMP s). Values higher than 0.8 are highlighted in boldface. Proximal Diaphysis Distal Global sf 0.718 0.879 0.660 0.844 AMP s 0.694 0.862 0.640 0.813 The inﬂuence of the mesh (C3D4 vs. C3D10) was analysed by comparing the correla- tion coefﬁcients in Table 6; in particular, the Spearman and the power Pearson coefﬁcients in the case GLT = 80 and using the constitutive model IsoH . The Pearson coefﬁcients 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 Y and the strain magnitudes (not shown), did not improve their correlations. The other types of ﬁt (linear and quadratic) also improved with C3D10, although to a lesser extent. It is noteworthy that the Spearman coefﬁcients were almost identical in both meshes. Mathematics 2022 ,10, 3367 20 of 29 Table 6. Inﬂuence of the FE mesh (C3D4 vs. C3D10) on the weighted Spearman and Pearson (power) correlation coefﬁcients obtained for GLT = 80 and constitutive model IsoH : for the maximum variable throughout the cycle, RM P, and for the amplitude throughout the cycle, RA P. The same predictors analysed for Spearman and Pearson coefﬁcients are compared here. Values higher than 0.9 (in Spearman) or 0.8 (in Pearson) are highlighted in boldface; negative values are in red. PredictorType of C3D4 C3D10 Correlation RM P RA P RM P RA P s1 Spearman 0.906 0.913 0.904 0.912 s3 Spearman 0.058 0.899 0.236 0.896 st Spearman 0.906 0.906 0.904 0.904 sc Spearman 0.895 0.895 0.893 0.894 sf Spearman 0.928 0.928 0.930 0.930 AMP s Spearman 0.922 0.921 0.923 0.922 svonMises Spearman 0.924 0.924 0.925 0.925 sTresca Spearman 0.923 0.923 0.925 0.924 so Spearman 0.742 0.913 0.734 0.911 #1 Spearman 0.431 0.429 0.499 0.496 #3 Spearman 0.431 0.443 0.484 0.495 AMP # Spearman 0.425 0.423 0.482 0.480 Y Spearman 0.722 0.721 0.709 0.709 Ytfor n = 2 Spearman 0.092 0.092 0.182 0.182 #o Spearman 0.254 0.464 0.287 0.499 AMP s Pearson (Power) 0.809 0.808 0.874 0.874 st Pearson (Power) 0.783 0.783 0.818 0.818 sc Pearson (Power) 0.692 0.693 0.741 0.742 sf Pearson (Power) 0.842 0.842 0.902 0.902 Y Pearson (Power) 0.006 0.006 0.011 0.011 Ytfor n = 1 Pearson (Power) 0.017 0.017 0.019 0.019 Ytfor 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 difﬁcult to distinguish. Instead, histograms are used to show the percentage of volume occupied by those elements whose values of AMP sor AMP #are within a given range. The elements of cortical ( r>1.2g/cm3) and trabecular bone ( r1.2g/cm3) 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 ( r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 of200–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 05101520253035% µ200–600600–10001000–1400 1400–1800 1800–2200>2200 0–200 AMP εcortical 200–6000102030405060708090 600–10001000–1400 1400–1800 1800–2200>2200% µε0–200 AMP σtrabecular 05101520253035404550 0–0.001 0.001–0.01 0.01–0.1 0.1–1 1–2 >2% MPa AMP σcortical 05101520253035404550 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 s) and strain (AMP #) for the elements of trabecular ( r1.2g/cm3) and cortical bone ( r<1.2g/cm3). The Y axis shows the percentage of the total volume occupied by the elements, whose AMP sor 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 coefﬁcients (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 coefﬁcients showed. The high correlation found between density and the von Mises, Tresca, absolute maximum principal stress (AMP s) and the ﬂuctuation of stresses ( sf) 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 r. 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é ( Yt) is even worse, as it yields negative correlation coefﬁcients. Mathematics 2022 ,10, 3367 22 of 29 The Spearman coefﬁcient is probably the most simple indicator that two variables are correlated by a monotonic function, since it does not assume any speciﬁc relationship as does the Pearson coefﬁcient. If R Sis 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 speciﬁc 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 speciﬁc variable should be further investigated. In this regard, the linear, quadratic or power functions tried for the Pearson coefﬁcients worked only moderately well. Only sfand AMP sshowed a value of R Pslightly over 0.8 (up to 0.861 for sfin the IsoJmodel). This means that R2>0.64, i.e., around 64% of the variation of the density, can be explained by sf(up to 74% in the IsoJmodel). Only some of the variables appearing in Table 1 were chosen for assessing its Pearson coefﬁcient, those having a high Spearman coefﬁcient plus SED and its related variables (Yt), for its common use as a mechanical stimulus in many models of bone remodelling. Obviously, if the Spearman coefﬁcients of SED and its related variables were low, the re- spective Pearson coefﬁcients could not improve them, but it is worth noting the very low values of R Pobtained 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 s) and one representing the strain-state (AMP #), was analysed by means of histograms (see Figure 6) that distinguish between cortical ( r>1.2g/cm3) and trabecular bone ( r1.2g/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 #2[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 #2[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 conﬁrmed by the high Spearman correlation coefﬁcients of most stress variables (see Table 1). The distributions and correlations of strain and stress variables would conﬁrm 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 Sin 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 ﬂuctuations of stresses throughout the cycle seems the more appropriate mechanical stimulus among those tested here, though the speciﬁc 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, Sr, is being forced by the Mathematics 2022 ,10, 3367 23 of 29 relationship between stiffness and density, E=F(r), in particular for those predictors derived from the stress tensor (for the sake of generality, we should write C=F(r), with Cbeing 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(r)) . It is well known from the literature that Fis a monotonically increasing function and given that the stiffer elements tend to withstand higher stress levels, we should expect a function such as s=H(E)to be monotonically increasing as well. Thus, if we write: s=H F(r) , (38) we can expect a monotonically increasing function relating sandr, in other words, a posi- tive correlation Srfor the predictors derived from the stress tensor. Therefore, it could seem that the function E=F(r)is forcing Srand it is so to some extent, especially for the Spearman correlation coefﬁcients. However, if E=F(r)were the only determinants of the correlation Sr, this correlation should depend on the constitutive model and it does not—signiﬁcantly, at least. Indeed, the increase in stiffness from IsoJtoIsoH 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 coefﬁcients or the AMP shistograms, 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 signiﬁcant effect on the correlations. However, there is another argument to support that Sris 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: s=H F(r),F , (39) where Fstands 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(s), thus leading to: S=G H F(r),F . (40) In view of Equation (40) the predictors, S, need not be predetermined only by the function F(r), not even the Spearman correlation coefﬁcients, as the function Gis not necessarily monotonic. This is conﬁrmed by the different Spearman coefﬁcients obtained for different predictors (see Tables 1 and 3). It should be noted that the same concepts that are behind the Gfunctions 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 Sr. We have obtained these correlations for every instant of the gait cycle separately. The stress and strain tensors obtained at each instant iwere used to assess Si e(recall Section 2.7) and these variables were correlated with reto give Ri P, a weighted Pearson correlation coefﬁcient for each instant i. The worst coefﬁcient between AMP sand den- sity throughout the cycle was Ri P= 0.48 (for the quadratic ﬁt and in the case GLT = 80, constitutive model IsoH and C3D4 elements), which is far from RM P= 0.804. The former is a very poor correlation. (Ri P)2= 0.23, i.e., only about 25% of the variance of density can be explained by AMP siand 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(r)is contributing to obtain the correlations Sr, at least for stress-based predictors, but only to some extent. There are other factors, not inﬂuenced 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 predict- ing 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 ﬁts (not shown) also improved with C3D10, but to a lesser extent. Notwithstanding, the C3D4 mesh had a sufﬁcient 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 coefﬁcients derived from both are also almost identical. For that reason, the use of C3D4 elements was justiﬁed, 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 inﬂuence of other activities can be expected. The joint loads were applied in a simpliﬁed way, as nodal-concentrated loads instead of loads distributed over the articular surface. This approximation is especially signiﬁcant in the knee reaction force and this could explain the lower correlation coefﬁcients 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 inﬂuence 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 conﬁrm 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 trabecu- lar 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 ﬂuctuation 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 Con- sejerí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íﬁco 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 ﬁnite 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. Conﬂicts of Interest: The authors declare no conﬂicts 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 s Absolute maximum principal stress MF Muscle force(s) HJF Hip joint force KJF Knee joint force Appendix A. Weighted Pearson and Spearman Correlation Coefﬁcients In the case of a linear correlation, the Pearson coefﬁcient between the bone apparent density rand the predictor varis weighted with the volume of the ﬁnite elements, which are the sample points. The weighted Pearson coefﬁcient is deﬁned as the weighted covariance divided by the weighted variances of both variables and determines the proportion of variance explained by the weighted linear ﬁt of the point cloud ( re,vare): RP=n å e=1[Ve(rer)(varevar)] sn å e=1[Ve(rer)2]n å e=1[Ve(varevar)2](A1) Mathematics 2022 ,10, 3367 26 of 29 with the weighted averages being , for example, in the case of the predictor var: var=n å e=1vareVe n å j=1Ve(A2) and with an analogous expresion for the weighted average of density, r.Veis the volume of the element e, used to weight the variable vare. The weighted Spearman correlation coefﬁcient, RS, is obtained through Equations (A1) and (A2) by simply replacing the pairs ( re,vare) with the pairs composed of their respective ranks, i.e., rank(re), rank (vare) . The weighted Pearson correlation coefﬁcient was also obtained for both a quadratic and a power ﬁt of the point cloud ( re,vare). In the ﬁrst case, the following function was ﬁtted to the point cloud in the least squares sense, after being weighted with the element volume, using the polyﬁtweighted library of MATLAB: ˜var(r) =Ar2+Br+C. (A3) In least squares ﬁtting the Pearson coefﬁcient is deﬁned as the quotient between the covariance of the ﬁtted and raw variables divided by the variances of both. In weighted ﬁttings Equation (A1) is modiﬁed as follows: RP=n å e=1[Ve(˜var(re)˜var)(varevar)] sn å e=1[Ve(˜var(re)˜var)2]n å e=1[Ve(varevar)2]. (A4) In the case of a power ﬁt, the function to be ﬁtted ( var=Krb) is transformed into a linear function by taking logarithms of both sides: logvar=K+blogr. (A5) Then, the Pearson coefﬁcient of a linear ﬁt (Equation (A1) ) is used by replacing the pairs ( re,vare) with (log re, log vare). Appendix B. Inverse Dynamic Problem Methodology The forces applied to the FE model were obtained solving an inverse dynamics prob- lem, which is brieﬂy described next (a more detailed description can be found for example in [38,39]). All the simulations were run in OpenSim. First, the inverse kinematics problem was solved. The input data in this problem are the trajectories of the markers, xext, obtained experimentally for the gait cycle by using a Vicon ®system. The output is the temporal evolution of the generalized coordinates, q. The inverse kinematics problem was solved using a least square pose estimator: minq å i2markerswikxexp ixi(q)k2! , (A6) where xi(q)is the position of the virtual marker i, which depends on the coordinates values, and wiis a marker weight taken from [ 17,21]. The results of the inverse kinematics problem were used as input to solve the inverse dynamics problem, by means of the classical equations of motion: M(q)¨ q+C(q,˙ q) +G(q) =ø, (A7) Mathematics 2022 ,10, 3367 27 of 29 where q,˙ qand ¨ qare the vectors of generalized positions, velocities and accelerations, respectively; Mis the mass matrix; Cis the vector of quadratic velocities; Gis the vector of external forces including gravitational forces and ground reaction forces and øis the vector of generalized forces regarding motor tasks. Finally, muscle forces ( Fmus) were estimated solving the following optimization problem: min J=n å m=1(am)p! s.t.n å m=1[amf(F0 m,lm,vm)]rm,j=tj, (A8) where nis the number of muscles in the model; amis the activation level of muscle m at a given time step; F0 mits maximum isometric force; lmits length; vmits shortening velocity; f(F0 m,lm,vm)its force-length-velocity relation; rm,jits moment arm about the jth joint axis and tjis the generalized force acting about the jth joint axis. The cost function J of Equation (A8) minimizes the sum of muscle activation squared (p = 2) as in [ 17]. As an example, Figure A1 shows the temporal evolution of the forces exerted by the iliacus and the lateral gastrocnemius obtained following this procedure. The temporal evolutions of all the muscles considered in the analysis are included in the Supplementary Materials. Figure A1. Temporal evolution of the components of the forces exerted by the iliacus and the lateral gastrocnemius throughout the gait cycle. AP: antero-posterior. CC: Cranio-caudal. LM: Lateral-medial. Finally, the bone-on-bone forces on the hip and the knee were calculated using the algorithm proposed by Steele [22]. Ri=Mi(y)¨ yi(åFmuscles +åFexternal +Ri1) (A9) In this equation, vector Ricontains the resultant forces and moments at joint ior bone-on-bone forces. The body distal to joint i,Bi, is treated as an independent body with known kinematics in a global reference frame. Thus, ¨ yirepresents the six dimensional vector of known angular and linear accelerations of Bi, while Mi(q)is the 66 mass matrix of Bi.Fexternal and Fmuscles represent the previously calculated forces and moments produced by external loads and musculotendon actuators, respectively. R(i1)represents the joint reaction load applied at the distal joint and was calculated in the previous recursive step. Our aim was to study bone-on-bone forces only at the hip and the knee joints, thus, Mathematics 2022 ,10, 3367 28 of 29 the generic name Riwill be replaced by HJF and KJF (acronyms for hip and knee joint forces, respectively) and expressed in the local frame attached to the femur. Figure A2 shows the temporal evolution of HJF and KJF throughout the gait cycle. These data are included in the Supplementary Material. Figure A2. Temporal evolution of the HJF and KJF throughout the gait cycle. AP: Antero-posterior. CC: Cranio-caudal. LM: Lateral-medial. BW: Body weight.