**Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A. Weighted Pearson and Spearman Correlation Coefficients**

In the case of a linear correlation, the Pearson coefficient between the bone apparent density *ρ* and the predictor *var* is weighted with the volume of the finite elements, which are the sample points. The weighted Pearson coefficient is defined as the weighted covariance divided by the weighted variances of both variables and determines the proportion of variance explained by the weighted linear fit of the point cloud (*ρe* , *vare*):

$$R\_P = \frac{\sum\_{\varepsilon=1}^{n} \left[ V\_{\varepsilon} \left( \rho\_{\varepsilon} - \overline{\rho} \right) \left( var\_{\varepsilon} - \overline{var} \right) \right]}{\sqrt{\sum\_{\varepsilon=1}^{n} \left[ V\_{\varepsilon} \left( \rho\_{\varepsilon} - \overline{\rho} \right)^{2} \right]} \cdot \sum\_{\varepsilon=1}^{n} \left[ V\_{\varepsilon} \left( var\_{\varepsilon} - \overline{var} \right)^{2} \right]} \tag{A1}$$

with the weighted averages being •, for example, in the case of the predictor *var*:

$$\overline{var} = \frac{\sum\_{\varepsilon=1}^{n} var\_{\varepsilon} \,\, V\_{\varepsilon}}{\sum\_{j=1}^{n} V\_{\varepsilon}} \tag{A2}$$

and with an analogous expresion for the weighted average of density, *ρ*. *Ve* is the volume of the element *e*, used to weight the variable *vare*.

The weighted Spearman correlation coefficient, *RS*, is obtained through Equations (A1) and (A2) by simply replacing the pairs (*ρe* , *vare*) with the pairs composed of their respective ranks, i.e., rank(*ρe*) , rank(*vare*) .

The weighted Pearson correlation coefficient was also obtained for both a quadratic and a power fit of the point cloud (*ρe* , *vare*). In the first case, the following function was fitted to the point cloud in the least squares sense, after being weighted with the element volume, using the *polyfitweighted* library of MATLAB:

$$
v \tilde{a}r(\rho) = A\rho^2 + B\rho + \mathbb{C}.\tag{A3}$$

In least squares fitting the Pearson coefficient is defined as the quotient between the covariance of the fitted and raw variables divided by the variances of both. In weighted fittings Equation (A1) is modified as follows:

$$R\_P = \frac{\sum\_{\varepsilon=1}^{n} \left[ V\_{\varepsilon} \left( v \overline{a} r (\rho\_{\varepsilon}) - \overline{v \overline{a} r} \right) (v a r\_{\varepsilon} - \overline{v a r}) \right]}{\sqrt{\sum\_{\varepsilon=1}^{n} \left[ V\_{\varepsilon} \left( v \overline{a} r (\rho\_{\varepsilon}) - \overline{v \overline{a} r} \right)^{2} \right]} \cdot \sum\_{\varepsilon=1}^{n} \left[ V\_{\varepsilon} \left( v a r\_{\varepsilon} - \overline{v a r} \right)^{2} \right]} \tag{A4}$$

In the case of a power fit, the function to be fitted (*var* = *K ρβ*) is transformed into a linear function by taking logarithms of both sides:

$$
\log var = K + \beta \log \rho.\tag{A5}
$$

Then, the Pearson coefficient of a linear fit (Equation (A1)) is used by replacing the pairs (*ρe* , *vare*) with (log *ρe* , log *vare*).

### **Appendix B. Inverse Dynamic Problem Methodology**

The forces applied to the FE model were obtained solving an inverse dynamics problem, which is briefly 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, **<sup>x</sup>ext**, 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:

$$\min\_{\mathbf{q}} \left( \sum\_{i \in markrs} w\_i \left|| \mathbf{x\_i^{exp}} - \mathbf{x\_i(q)} \right||^2 \right),\tag{A6}$$

where **xi**(**q**) is the position of the virtual marker *i*, which depends on the coordinates values, and *wi* is 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:

$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}}) + \mathbf{G}(\mathbf{q}) = \mathbf{o},\tag{A7}$$

where **q**, **q˙** and **q¨** are the vectors of generalized positions, velocities and accelerations, respectively; **M** is the mass matrix; **C** is the vector of quadratic velocities; **G** is 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:

$$\begin{aligned} \min \left( \left. f = \sum\_{m=1}^{n} (a\_m)^p \right| \right) \\ \text{s.t.} \quad \sum\_{m=1}^{n} \left[ a\_m \cdot f(F\_m^0, l\_m, \upsilon\_m) \right] \cdot r\_{m,j} = \tau\_j \end{aligned} \tag{A8}$$

where *n* is the number of muscles in the model; *am* is the activation level of muscle *m* at a given time step; *F*0*m* its maximum isometric force; *lm* its length; *vm* its shortening velocity; *f*(*F*0*m*, *lm*, *vm*) its force-length-velocity relation; *rm*,*j* its moment arm about the *j*th joint axis and *τj* is the generalized force acting about the *j*th 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].

$$\mathbf{R\_i} = \mathbf{M\_i(y)}\vec{\mathbf{y\_i}} - \left(\sum \mathbf{F\_{muscle}} + \sum \mathbf{F\_{external}} + \mathbf{R\_{i-1}}\right) \tag{A9}$$

In this equation, vector **Ri** contains the resultant forces and moments at joint *i* or 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, **y¨i** represents the six dimensional vector of known angular and linear accelerations of *Bi*, while **Mi**(**q**) is the 6 × 6 mass matrix of *Bi*. **Fexternal** and **Fmuscles** represent the previously calculated forces and moments produced by external loads and musculotendon actuators, respectively. **<sup>R</sup>**(**<sup>i</sup>**−**<sup>1</sup>**) 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, the generic name **Ri** will 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 Materials.

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