Using a Truss-Inspired Model with the Uniform Strength Optimization Theory to Predict Spongy Bone Geometry in Proximal Femur

This paper presents a new naïve approach for simulating bone remodeling process. It is based on the uniform strength theory of optimization and employs a truss-like model for bone. The truss was subjected to external loads including 5 point loads simulating the hip joint contact forces and 3 muscular forces at the attachment sites of the muscles to the bone and the rest are reactions of ligaments. The strain in the links was calculated and the links with high strains were identified. The initial truss is modified by introducing new links wherever the strain exceeds a prescribed or critical value. The critical value was assumed to be equal to an average of the absolute value of strains in the initial model. Each link which undergoes a high strain is replaced by several new links by adding new nodes around it using Delaunay method. Introducing the new links to the truss, which is conducted according to a weighted arithmetic mean formula, will strengthen the structure and reduce the strain within the respective zone. This procedure was repeated for several times. Convergence was achieved when there were no critical links remaining. This method was used to study the 2D shape of proximal femur in the frontal plane and provided results which are in a fairly good agreement with CT image of the human proximal femur.


INTRODUCTION
An early hypothesis about the dependence of the structure and form of bones on their mechanical function was proposed by Galileo in 1638 [1]. The nature of this dependence was first described in a semiquantitative manner by Wolff [2], who stated that every change in function of a living bone is followed by adaptive changes in its internal architecture and shape. Several mechanisms have been proposed to relate changes in mechanical loads to the adaptive responses in bone, including: piezoelectric and streaming potentials [3,4], mechanical fatigue microdamage [5][6][7][8][9][10][11], and extra-cellular fluid pressure gradient effects on bone cells [12,13]. Experimental evidence can be found in support of each of the above-mentioned mechanisms. It is generally accepted that bone growth, maintenance, and degeneration are biochemically regulated processes that are influenced by mechanical functions [14,15]. Remodeling theories can be categorized into phenomenological, mechanistic, and optimization models [16]. Generally speaking, a remodeling theory tries to find a relation between bone apparent density and the mechanical stimuli applied on the bone. Also, a remodeling theory can aim at finding a relation between the mechanical stimuli and the bone micro-structural pattern, i.e. the fabric tensor [17]. Many theories, mostly phenomenological, have been proposed to explain the process of bone remodeling [19][20][21][22][23][24][25][26][27][28]. For example, the apparent density has been used to characterize the internal morphology of the bone and the strain energy density assumed as the stimulus for bone adaptation [19]. These methods, however, are often based on complicated mathematical formulations and/or require extreme computation costs. A high-order nonlinear equation of bone remodeling was utilized in [21] and the influence of nonlinearity on adaptation was studied. Prediction of external shapes and internal density distributions in the proximal femur was also investigated using these non-linear equations [22].
It is expected that the application of different boundary conditions to the Finite Element model (FEM) will have an effect on the density distribution of bone. Recent studies have shown significant differences in the femoral strain distribution under varying boundary conditions [29][30][31][32]. So, it is expected that different boundary conditions will have a direct effect on the density distribution in the proximal femur. It is not clear, however, how these different loading configurations can affect the outcome of the remodeling simulation. It has been also argued [19,33] that representation of all muscle forces about a joint is not needed.
The goals of this study are to find the external shape of proximal femur and its density distribution using a novel, computationally efficient remodeling method. The effects of different boundary conditions on the density distribution will also be investigated.

ANALYTICAL METHOD
In order to find the external shape of bone and the bone density distribution, a topology optimization method was employed. The initial design domain is a uniform rectangular plate which occupies a larger space than the anticipated final shape of the bone structure. The rectangle is meshed by 156 link elements ( Fig. 1) and subjected to the actual loadings found on the proximal femur ( Fig. 2 and table I). This consists of muscle contributions from the gluteus medius, gluteus minimus, piriformis, and hip joint contact force [22 and 29].
Strain in each link is calculated using an algorithm written in Matlab. A critical strain threshold is established and links with strains above this critical value ( c ) are identified. Additional nodes are added near any critical links, i.e. the links which have strains above c , in order to stiffen the structure in the vicinity of the high strain links.
The strain is then recalculated in the structure. At each step critical links are identified and additional nodes are added. Convergence is a function of the critical strain; thus, a judicious choice of critical strain is needed to minimize computational time. We have found that the average of the absolute value of strains in the initial model is a good value for the critical strain threshold. In fact, this is the advantage of the present method in which strain tensor of each region is reduced to a scalar. All the strains of nodes together make a strain vector which is here a linear combination of displacement vector of the whole structure. Since, the essential, i.e. Dirichlet or displacement, boundary condition of the problem is homogeneous; the displacement vector is inversely proportional to the Young's modulus of the links which has been considered identical throughout the structures. This can be seen in Eq. (1): 22 21 12 11 (1a) , 1 1 11 In Eq. (1a), 0 is a null sub-vector assigned to the supported nodes and  Since, the stiffness matrix is linearly proportional to the Young's modulus (E), then Eq. (1b) implies that all displacement values are inversely proportional to Young's modulus (tilde symbols are independent of E). Consequently, the strain vector which is a linear combination of displacement vector is also inversely proportional to the Young's modulus. Stress vector, on the other hand, is Young's modulus times strain vector, resulting in independence of stress vector to Young's modulus. The above statements have been restated in terms of mathematical symbols as follows: There exist two optimization theories: uniform strength theory and trajectorial architecture theory [34]. In the first approach, structure is made in a way that under loading condition of interest, the maximum allowable stress all over the structure is identical. In the latter approach, however, the structure is made such that material appears only in the trajectories of loading transfer in a particular loading condition. In this study, the first approach has been considered. As already stated, the stress vector is independent of Young's modulus. Therefore, to achieve an optimized structure, we do not need to have bone's Young's modulus in this method. For any arbitrary value of Young's modulus, a vector of strains and then stresses are obtained. Uniformity of stress distribution necessitates the uniform distribution of strain. The obtained strains span from a minimum value, min to a maximum value, max . By trial and error, it is concluded that if we use the average of min and max then a convergence to optimal shape has an acceptable rate. If one chooses a smaller amount, then more iteration might be needed to achieve the optimum shape. Also, if a greater amount is chosen for the strain, then more regions would undergo strengthening and the structure would become more elaborated. The elaborated structure means one with too many links. This, in turn, implies that in order to obtain displacements in the latter steps a too big system of equations should be solved and, consequently, more overall computational cost is spent in addition to oversensitivity to the loads. The over-sensitivity, on the other hand, is because much bigger region is influenced by each particular load. In other words, when a particular load is applied, then more proportion of structure is considered as critical and hence subject to reinforcement. It is clear that this process contradicts with the aim of optimization.
New nodes are added by means of a weighted method in which the added nodes are closer to the critical link. The nearer the new nodes to the weak link, the stronger link will result. Not to influence the uniformity of nodes distribution, we have used the following weighted scheme. Consider that the link 1-2 in Fig. 3 is a critical link, and two nodes are going to be added to it. The location of the new nodes, i.e. nodes 5 and 6 in Fig. 3, are found by using the following relations:  The new nodes are attached to the initial nodes with new links using the Delaunay method. The Delaunay triangulation is one of the most popular methods for generation of unstructured meshes. It is composed of two phases: placement of the mesh vertices, and triangulation. The Delaunay algorithm determines non-overlapping triangles that fill the convex hull of the input points such that every edge is shared by at most two triangles, and the circumcircle of every triangle contains no other input points [35]. In this way, a sub-structure will be resulted in which stiffness between nodes 1 and 2 has been enhanced. For example, when two triangles 123 and 124 in Fig. 3 are equilateral, the strengthened sub-structure obtained by this method between nodes 1 and 2 is 1.48 times stiffer than initial sub-structure. In other words, if the right structure in Fig. 3 is supported in node 2 and is stretched in node 1, the displacement of node 1 will be 1.48 times less than displacement of node 1 in the left structure of Fig. 3. Figure 4 shows the model after convergence. The number of nodes after convergence is approximately 1500 and this means that the code will not add any new nodes. As a final step, a "density matrix" is developed to find the density distribution in the model. The model is divided into 100×100 small rectangles. The density is found by dividing the mass of the trusses in the rectangle by its volume, which is assumed to have unit thickness. The mass of the trusses in each rectangle is found by multiplying the volume of trusses by 1/ max =1.0gr/cm 3 which is a typical bone density [23]. By assigning white color to the high density areas and dark color to the low density areas, a gray-scale image is generated which is shown can be seen in Fig. 5. Figure 6 is a C-T image of human proximal femur. The white regions correspond to high density bone and black regions to either low bone density or bone porosity. An interesting point about Fig. 5 is that the highest density can be seen in zone 1 of our model which corresponds to the cortical shell of the proximal femur (see Fig. 6). A fairly good agreement between Figs. 5 and 6 can be seen.

RESULTS AND DISCUSSIONS
The convergence rate in this model is very high, with an appropriate value for critical strain. As it is shown in Fig. 7, the convergence is achieved after only 10 runs.  Table I) in the remodeling procedure, different density distributions are achieved. For example, by decreasing or increasing the hip joint contact forces by 10% (P1, P2, P3, P4, and P5), the density distribution will be altered (Figs. 8 and 9). The average change in the density is about -0.1% and 2%, respectively. The reason why changes in average density as a result of increasing and decreasing of loads is not equal, or at least comparable, is that the effect of bone resorption has not been taken into account in this study, and only bone formation process has been considered in our model. Results of these simulations put an emphasis on the well accepted point in the biomechanics of bone that magnitude and also direction of loads applied on bone (from muscles, tendons, or ligaments) has a significant influence on the bone apparent density, and also its microstructure.

CONCLUSIONS
The whole topology of a bone structure is made up of its external shape and internal density distribution. In the present study, we put forward and fulfill the prediction of the external shape of the human proximal femur, and its density distribution by a fast and simple method. In this study, unlike most of the existing models of bone remodeling which are on a continuum mechanics ground, a discontinuous, truss like model to analyze the bone remodeling process in the human proximal femur is used. The philosophy behind our model is the very well accepted point of correspondence between the level of external loads and the bone strength and density. Uniform strength optimization method was adopted in this research which resulted in a fairly fast convergence time (about 20 minutes). Using our method, it is found that the density distribution would be fairly similar to the real human proximal femur (Figs. 5 and 6). Furthermore, it was shown that by decreasing the external loads on the femur, there will be a reduction in the bone density too (Fig. 8).
Considering that the real proximal femur is a three dimensional structure, a 3-D modeling of the bone with the same method is intended for our future work. Moreover, in order to increase the clarity of some of our results, e.g. Figs.5 and 8, more research is in progress.